Atomic characteristics and high-order harmonic spectra:
extension of ab-initio numerical calculations to larger
systemsby
Ran Brynn Reiff
B.A., University of Colorado Boulder, 2015
M.S., University of Colorado, 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
2021
Committee Members:
Andreas Becker, Chair
Agnieszka Jaroń-Becker
John Bohn
Bengt Fornberg
Margaret Murnane

Reiff, Ran Brynn (Ph.D., Physics)
Atomic characteristics and high-order harmonic spectra: extension of ab-initio numerical calculations to larger systems
Thesis directed by Prof. Andreas Becker

High-order harmonic generation (HHG) is a highly nonlinear process where an electron driven
by an intense, typically infrared, laser pulse ionizes, accelerates in the laser field and recombineswith the parent ion, emitting a photon with frequency many times that of the driving field. The
coherent radiation emitted from many atoms in the focus of the driving laser leads to the generationof an ultrashort ultraviolet or x-ray laser pulse. Numeric solution of the time-dependent Schrödinger
equation (TDSE) accurately describes the single-electron process, including the low-energy regionof the emission spectrum for which the atomic potential and excited states play a significant role.
The coherent macroscopic response (e.g., from a gas jet) involves the emission of a very large
number of atomic radiators (e.g., 1015 ) that each interacts with the driving laser at different peakintensities and carrier-envelope phases. Full ab-initio simulation of the macroscopic response using
exact numeric single-electron calculations are not feasible.
In this thesis, we present results of three projects related to high-order harmonic generation with a focus on the below- and near-threshold regime, in which excited states of the atomand the specific form of the atomic potential play an important role. First, we present a reproducible ab-initio method to produce benchmark tests between calculations based on the timedependent Schrödinger equation (TDSE) in the single-active-electron approximation (SAE) andtime-dependent density functional theory (TDDFT) in the highly nonlinear multiphoton and tunneling regime of strong-field physics. As key to the benchmark comparison we obtain an analytic
form of SAE potentials based on density functional theory, which we have applied for different
atoms, ions, and molecules. Using these potentials, we find remarkable agreement between the
results of the two independent numerical approaches (TDDFT and SAE-TDSE) for the high-orderiii
harmonic yields in helium, demonstrating the accuracy of the SAE potentials as well as the predictive power of SAE-TDSE and TDDFT calculations for the nonperturbative and highly nonlinear
strong-field process of high harmonic generation in the ultraviolet and visible wavelength regime.
Next, we investigate resonance enhancement of near-threshold HHG, identifying similar resonance effects in broad parameter regimes, including hydrogen driven by near- and mid-infrared
pulses as well as helium with 400 nm lasers. We differentiate behavior that may be explicablethrough semi-classical trajectory models (generally at higher intensity) from features which are notconsistent with these models.
In the last part of the thesis, we develop a macroscopic description of high-order harmonic
radiation resulting from the interaction of atomic systems with an intense laser pulse using abinitio solutions of the time-dependent Schrödinger equation (TDSE). We show that for this highly
nonlinear process, interpolation can be performed across laser intensity for a given wavelength,
limiting the number of full time-dependent Schrödinger equation calculations to about one hundred.
The significantly reduced computational time as compared to more sophisticated methods opens apath toward the extension of macroscopic high harmonic calculations based on ab-initio microscopic
results to more complex targets and interactions. We investigate the near-threshold regime of thespectra, showing that the degree of coherence of the off-harmonic radiation generated during thepulse is much lower than that of the harmonics, but this radiation extends to larger divergenceangles than the harmonic signals – providing the option for separation of the different signals inthis part of the spectrum. Finally, we analyze high harmonic generation from the spatial phase
distribution for broadband Gaussian pulses with a negative Porras factor, showing an interference
pattern in the angular distribution of below- and near-threshold harmonics, which is not present
for the monochromatic Gouy phase distribution.

Dedication

To my wife, Lorien.

Acknowledgements

First, my research mentors and colleagues. The Ultrafast AMO Theory group has welcomedme since I began an undergraduate research project in 2014. Profs. Andreas Becker and Agnieszka
Jaroń-Becker introduced me to research, and have since taught me countless things about physics,
research, academia, and teaching. Senior grad students Dr. Michelle Miller and Dr. Cory Goldsmith were particularly important to my onboarding, both scientific and social. Later, I was lucky
enough to grow and work alongside Erez Shani, Dr. Joel Venzke, and Dr. Tennesse Joyce. We
complemented each other well with our different perspectives and approaches, leading to many
productive debates at the white board. Newer graduate students (Spencer Walker, Yonas Gebre,
Andres Mora, and Bejan Ghomashi) and undergraduates (Lucas Kolanz and Zetong Xue) challenged me to rethink and repackage my knowledge, expanding it in the process. As a postdoc at
JILA and later, Dr. Carlos Hernandez-Garcia provided crucial expertise and I am indebted to his
aid on my work understanding macroscopic HHG.
Teachers beyond physics have played a vital role in cultivating ways of thinking and communicating necessary for research. I do not have the words (or space) to convey the impact each of these
individuals had on my life and worldview: Peg Bemis (literature), Marco de Martino & Anissa Butler (history, sociology, and writing), Kira Hall (sociolinguistics), Christian Kopff (classics), Valerie
Otero (pedagogy and systemic change).
My research has been bolstered by the community and ways of thinking that I developed
doing equity work. Without that support and motivation, I would not have completed this thesis.
I am especially grateful to Dr. Miller for bringing me into this work via CU-Prime, from which I

viextended my involvement into many other groups (R-Cubed, JEDI, APS-IDEA, and GASP). The
friendship and shared mission of improving physics have continually buoyed me when I struggled
in research and life.
Finally, I must thank my mother, my wife, and my friends, who encouraged the persistent spirit of lifelong learning, thoughtful collaboration, and critical thinking necessary for research.

Contents

Chapter
1 Introduction
1.1

1.2

1.3

1.4

1.5

1

Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.1.1

Single- and multi-photon absorption . . . . . . . . . . . . . . . . . . . . . . .

1

1.1.2

Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.1.3

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

2Strong field ultrafast physics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.2.1

Laser pulses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.2.2

Tunnel and multiphoton ionization . . . . . . . . . . . . . . . . . . . . . . . .

5

AC Stark shift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.3.1

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

1.3.2

Free-induction decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

Numeric methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4.1

Semi-classical trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

1.4.2

Finite difference and eigenstate problems . . . . . . . . . . . . . . . . . . . . 15

1.4.3

Time-dependent Schrödinger equation . . . . . . . . . . . . . . . . . . . . . . 16

Tools for analyzing HHG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.5.1

Temporal windowing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

1.5.2

Spectral filtering and temporal profiles . . . . . . . . . . . . . . . . . . . . . . 21

1.5.3

Wavelet analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

viii
1.5.4
1.6

Parameter scans . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

Outline of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2 Single-Active-Electron Potentials
2.1

2.2

2.3

Constructing SAE potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.1.1

Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.1.2

Notes on implementing DFT . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.1.3

Analytic fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

2.1.4

Ionization and excitation energies . . . . . . . . . . . . . . . . . . . . . . . . . 38

HHG with SAE potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.2.1

SAE-TDSE and TDDFT methods . . . . . . . . . . . . . . . . . . . . . . . . 40

2.2.2

Comparison of SAE-TDSE and TDDFT results . . . . . . . . . . . . . . . . . 41

2.2.3

Multiple SAE-TDSE for larger atoms . . . . . . . . . . . . . . . . . . . . . . 47

2.2.4

Notes on SAE for large atoms . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

2.2.5

Effects of manual SAE parameter modification . . . . . . . . . . . . . . . . . 49

Molecular SAE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.3.1

2.4

27

Potentials for other molecules . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.4.1

Usage in other work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

3 Resonance-Enhanced High-order Harmonic Generation

58

3.1

Spectral features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

3.2

Resonantly-initiated quantum trajectories . . . . . . . . . . . . . . . . . . . . . . . . 63

3.3

Identifying resonances with temporal profiles . . . . . . . . . . . . . . . . . . . . . . 65

3.4

Radiation at other driving intensities . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

3.5

CEP and the trajectory model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.5.1

3.6

Varying pulse envelope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

Deviations from predicted resonance times . . . . . . . . . . . . . . . . . . . . . . . . 87

ix
3.6.1

Deviations in radiation features around 9-photon resonance to 3p state . . . . 89

3.7

Short and long pulses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

3.8

Appearance in other parameter regimes . . . . . . . . . . . . . . . . . . . . . . . . . 101

3.9

3.8.1

400 nm helium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

3.8.2

1600 nm hydrogen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
3.9.1

Future prospects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

4 Interpolation and Macroscopic Response

112

4.1

Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

4.2

Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

4.3

4.2.1

Interpolation across intensity . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

4.2.2

Variation of CEP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

Macroscopic summation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
4.3.1

4.4

Test case: Lattice grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128

Atomic hydrogen gas jet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
4.4.1

Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

4.4.2

Near-threshold harmonics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

4.4.3

Appearance of resonance features in macroscopic response . . . . . . . . . . . 139

4.5

Effects of macroscopic phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

4.6

Other investigations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150

4.7

4.6.1

Hydrogen at 1600 nm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

4.6.2

Helium at 400 nm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159

5 Summary and Perspectives

164

Bibliography

167

List of Tables

Table
1.1

Description of Hartree atomic units, including reference constants and conversion tocommon units. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.1

3

Fitting results for charge parameters C0 and Zc in SAE potentials for rare gas atoms
helium, neon, and argon. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

2.2

Parameters for VSAE (Eq. (2.14)) obtained in analytical fits of single-active-electron
potentials for various atoms and ions. (The upper part of this Table also appears in[1].) (*) For Rb, Sr, Cd, and Xe, see Sec. 2.2.4.

2.3. . . . . . . . . . . . . . . . . . . . 36

Comparison of first ionization energies: Experimental data [2], DFT calculations,and analytical potential fit (SAE) (Taken from [1]).

2.4. . . . . . . . . . . . . . . . . . 38

Same as Table 2.3 but for ionization energies of inner orbitals in neon (Taken from[1]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

2.5

Comparison of experimental data [2] and predictions using the analytical potential
fit (SAE) for ionization energies of excited states in helium (Taken from [1]). . . . . 39

2.6◦ )X argon excited states (Taken from [1]). . . . . 39
Same as Table 2.5 but for 3p5 (2 P1/2

xi
2.7

Modifications of helium SAE potential keeping IP = 0.948. In the first row, potentials are labeled by which terms are included (Coulomb, Yukawa, Exponential)
and which term is modified from the base parameters. The second row defines theform of the modified potential, with specific parameter values given in the following
rows. Finally, the eigenstate energies are given through E4f . The potential withouta Coulomb term (YE-a) does not support bound states above |1si. We note theconsistency of state energies (and, as will be seen, HHG) for the first three potentials. 51

2.8

Parameters for molecular VSAE (Eq. (2.24-2.26)) obtained in analytical fits of singleactive-electron potentials for various molecules.

3.1. . . . . . . . . . . . . . . . . . . . 56

Identified resonances, delays, and line styles. We note that due to l-independent
energies in hydrogen atoms, some of the resonances can be ascribed to different
states with the same l parity, namely: 3s and 3d, 4s and 4d, 6p and 6f or 6h. . . . . 74

3.2

Classification of two-parameter scans keeping 9ω → 3p resonance time Tr fixed. For
each scan, two parameters are fixed at their base values and two are varied: oneparameter is given as a range of equally-spaced values, while the other is determinedby solving Eq. (3.3) analytically or numerically. The final four columns describe therange and trends of dimensionless parameters across the scans. . . . . . . . . . . . . 90

4.1

Comparison of parameters for hydrogen driven by near-IR and mid-IR lasers andfor helium driven by a 400 nm pulse. A single example intensity is chosen from therange of each set. We note that the similarity of dimensionless parameters for the
800 and 400 nm cases was by design.

. . . . . . . . . . . . . . . . . . . . . . . . . . 151

List of Figures

Figure1.1

Diagram of ionization processes. (black) Field-free atomic potential and 1s boundstate energy. (blue) Tunnel ionization: Potential distorted by external field, resulting
in tunnel exit x0 . (purple) Three-photon ionization resulting in excess kinetic energy
T0 .

1.2

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

Diagram of the AC Stark shift as a function of laser intensity for a wavelength of
800 nm. (solid lines) Excitation and (dashed line) ionization energies are given inunits of ω. As the intensity increases, N -photon ionization channels close (e.g., the9-photon channel at intensities above about 6 × 1012 W/cm2 ) and excited states shiftacross resonances.

1.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9Diagram of the time-dependent AC Stark shift. Near the peak of the pulse, theeffective ionization potential (top dashed line) rises above the energy of five photons(red arrows), while the 2s state energy (top solid line) shifts briefly into 4-photon
resonance.

1.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9Excited state distribution as function of quantum numbers n (vertical axis) and l(horizontal axis) at the end of 800 nm, 20 cycle pulses with sine-squared envelope
and peak intensities: (a) I0 = 3.4 × 1013 W/cm2 , (b) I0 = 6.0 × 1013 W/cm2 , (c)I0 = 8.6 × 1013 W/cm2 , and (d) I0 = 11.2 × 1013 W/cm2 . Left (right) column
corresponds to cases in which the Rydberg states are resonant with an even (odd)number of photons. (Taken from [3].). . . . . . . . . . . . . . . . . . . . . . . . . . 11

xiii
1.5

Diagram of a high-order harmonic generation spectrum marked by a plateau of radiation at similar intensities.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.6

Diagram of the 3-step model of high-order harmonic generation. (Taken from [4].) . 13

1.7

Radiation spectra generated during the pulse (red dashed lines) and those including
line emissions after the pulse (black solid lines) are shown for hydrogen driven by a
20 o.c. sine-squared, 800-nm pulse with peak intensities: (a) I0 = 3.4 × 1013 W/cm2 ,(b) I0 = 6.0 × 1013 W/cm2 , and (c) I0 = 16.4 × 1013 W/cm2 . The vertical gray lines
show the field free energy differences between the np energy levels (up to 14p) andthe 1s ground state. (Taken from [3].) . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.1

Results of DFT calculations for (a) electron density and (b) exchange correlation
potential of the helium atom (Taken from [1]). . . . . . . . . . . . . . . . . . . . . . 35

2.2

Same as Fig. 2.1 but for the argon atom (Taken from [1]).

. . . . . . . . . . . . . . 35

2.3

Relative error between results of the analytical fit and DFT calculations for the
single-active-electron potentials for (a) helium and (b) argon (Taken from [1]).

2.4. . . 37

Comparison of TDDFT results (dashed red line) and SAE-TDSE (solid black line)
results for HHG spectra generated in the interaction of a helium atom with a 4-cycle
laser pulse at peak intensity of 1015 W/cm2 and central wavelengths of (a) 267 nm,(b) 400 nm, and (c) 600 nm. Vertical dashed lines denote the ionization potential
IP and expected HHG cut-off energy for each calculation. (Taken from [1].)2.5. . . . 42

Results of SAE calculations for high-order harmonic spectra generated from (a) 2pshell and (b) 2s-shell in neon by a 4-cycle laser pulse at peak intensity of 1014 W/cm2
and wavelength of 800 nm. The vertical dashed line indicates the 2s − 2p transitionenergy at which an artificial peak occurs in the spectra from the two sub-shells,
which is removed once the the amplitudes for both processes are added (Panel (c)).
(Taken from [1]).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

xiv
2.6

Results of SAE calculations for high-order harmonic spectra generated in the interaction of a helium atom with a 4-cycle laser pulse at peak intensity of 1015 W/cm2
and central wavelength of 800 nm. Each line represents results for an SAE potential
defined in Table 2.7. The black line corresponds to the original potential as well as
two other potentials with indistinguishable spectra. Vertical dashed lines denote theshared ionization potential (IP = 0.948 a.u.) and HHG cut-off (8.0 a.u.). Results
have been normalized to the intensity at the cut-off. . . . . . . . . . . . . . . . . . . 52

2.7

Top row: Comparison of the analytical SAE potential fit (dashed red line) and thenumeric DFT potential (solid black line) (a) along the molecular axis (ρ = 0) and(b) for a line bisecting the molecular axis at the origin (o = 0). Bottom row: (c)Relative error between numeric DFT potential and analytic SAE fit (Eq. (2.24)).
Note that the error, within a range of ≈ 1.5 a.u. from the nuclei, is ≤ 5% (Taken
from [5]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

3.1

High harmonic spectrum generated in atomic hydrogen, exhibiting a rich structurein the near-threshold region of the spectrum between harmonics 5 and 12. The
hydrogen ionization threshold of 0.5 a.u. = 8.77ω, where ω = 0.057 a.u. corresponds
to the central driver wavelength of 800 nm. Other laser parameters: sin2 pulse of 20
optical cycles with peak intensity of 3.8 × 1013 W/cm2 . Note that the intensity ofthe spectrum is scaled by the signal at the fundamental frequency. . . . . . . . . . . 60

3.2

HHG spectra as a function of peak intensity for (top) radiation during the pulse (asin Fig. 3.1) and (bottom) for propagation of the wavefunction over a period of fourtimes the pulse duration after the end of the pulse. Results for each intensity havebeen normalized to the maximum signal in the region between 8.75ω and 9.5ω. Red
dashed lines denote field-free excited state transition energies.

. . . . . . . . . . . . 61

xv
3.3

HHG spectra as (top) a function of the driving wavelength at Nτ = 20 o.c., and(bottom) a function of pulse duration at λ = 800 nm. For both plots, driving intensity I0 = 3.8 × 1013 W/cm2 . Results include spectral contributions (free induction
decay) for propagation of the wavefunction over a period of four times the pulse duration after the end of the pulse. Results for each wavelength (top) or each duration(bottom) have been normalized to the maximum signal in the region between 8.75ω
and 9.5ω. Red dashed lines denote field-free excited state transition energies. . . . . 62

3.4

HHG spectral (top) intensity and (bottom) phase as a function of peak intensity forradiation during the pulse (other parameters as in Fig. 3.1).

3.5. . . . . . . . . . . . . 64

Temporal profile of radiation within 1ω of the 9th harmonic as a function of pulse
duration in atomic hydrogen for an 800 nm, 3.8 × 1013 W/cm2 pulse. For each
duration the results have been normalized to a maximum of 1 for sake of comparison.
A vertical dashed line denotes the peak of the laser and a horizontal dashed line
identifies the pulse parameters used as a base in this Chapter (i.e., 20-cycle duration). 68

3.6

Same as Fig. 3.5, but overlaid with lines denoting: (top) resonance times for 7(magenta) and 8- or 9- (white) photon resonances (there are no resonances with asmaller or larger number of photons involved for these parameters) for 2s − 6p, and(bottom) 9ω → 3p resonances delayed by 0.5 (solid), 1.8 (dashed), and 3.0 (dotted)
optical cycles.

3.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

Temporal profile of radiation within 1ω of the 9th harmonic as a function of driving
pulse intensity in atomic hydrogen for an 800 nm, 20 o.c. pulse. For each intensity the
results have been normalized to a maximum of 1 for sake of comparison. A verticaldashed line denotes the peak of the laser and a horizontal dashed line identifies thepulse parameters used as a base in this Chapter (i.e., 3.8 × 1013 W/cm2 intensity).
Panels show identical radiation data, but the lower is overlaid with lines denoting
9ω → 3p resonances delayed by 0.5 (solid), 1.8 (dashed), and 3.0 (dotted) optical
cycles.

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

xvi
3.8

Temporal profile of radiation within 1ω of the 9th harmonic as a function of driving
wavelength in atomic hydrogen for a 20 o.c. pulse at 3.8 × 1013 W/cm2 . For eachwavelength the results have been normalized to a maximum of 1 for sake of comparison. A vertical dashed line denotes the peak of the laser and a horizontal dashedline identifies the pulse parameters used as a base in this Chapter (i.e., 800 nm wavelength). Panels show identical radiation data, but the lower is overlaid with lines
denoting 9ω → 3p resonances delayed by 0.5 (solid), 1.8 (dashed), and 3.0 (dotted)
optical cycles.

3.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Temporal profile of radiation within 1ω of the 9th harmonic as a function of atomic
nuclear charge Z in a hydrogen-like atom for an 800 nm, 20 o.c. pulse at 3.8 × 1013
W/cm2 . For each Z the results have been normalized to a maximum of 1 for sakeof comparison. A vertical dashed line denotes the peak of the laser and a horizontal
dashed line identifies the pulse parameters used as a base in this Chapter (hydrogen,i.e., Z = 1). Panels show identical radiation data, but the lower is overlaid with lines
denoting 9ω → 3p resonances delayed by 0.5 (solid), 1.8 (dashed), and 3.0 (dotted)
optical cycles.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

3.10 HHG spectra as a function of peak intensity for radiation during the pulse. Resultsfor each intensity have been normalized to the maximum signal in the region between
8.75ω and 9.5ω. Red dashed lines denote field-free excited state transition energies.

75

3.11 Temporal profile of radiation within 1ω of the 9th harmonic as a function of pulse
duration in atomic hydrogen for an 800 nm pulse at (top) 2.0 × 1013 W/cm2 and(bottom) 3.0 × 1013 W/cm2 . For each duration the results have been normalized toa maximum of 1 for sake of comparison. A vertical dashed line denotes the peak of
the laser and a horizontal dashed line identifies the pulse duration used as a basein this Chapter (i.e., 20 o.c.). Lines are drawn for delayed resonances as outlined inTable 3.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.12 Same as Fig. 3.11, but for (top) 5.6 × 1013 W/cm2 and (bottom) 7.0 × 1013 W/cm2 .

79

xvii
3.13 Temporal profile of radiation within 1ω of the 9th harmonic as a function of driving
pulse intensity in atomic hydrogen for an 800 nm, 20 o.c. pulse. For each intensity the
results have been normalized to a maximum of 1 for sake of comparison. A verticaldashed line denotes the peak of the laser and horizontal dashed lines identify theintensities selected above. (lower plot) Lines are drawn for delayed resonances as
outlined in Table 3.1.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

3.14 Temporal profile of radiation within 1ω of the 9th harmonic as a function of driving
pulse wavelength in atomic hydrogen for a 20 o.c. pulse at (top) 2.0×1013 W/cm2 and(bottom) 3.0 × 1013 W/cm2 . For each wavelength the results have been normalized
to a maximum of 1 for sake of comparison. A vertical dashed line denotes the peakof the laser and a horizontal dashed line identifies the laser wavelength used as a basein this Chapter (i.e., 800 nm). Lines are drawn for delayed resonances as outlined inTable 3.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.15 Same as Fig. 3.14, but for (top) 5.6 × 1013 W/cm2 and (bottom) 7.0 × 1013 W/cm2 .

83

3.16 Temporal profile of radiation within 1ω of the 9th harmonic as a function of driving
pulse CEP in atomic hydrogen for an 800 nm, 20 o.c. pulse at (top) 3.8×1013 W/cm2 ,and (bottom) 5.6 × 1013 W/cm2 . A vertical dashed line denotes the peak of the laserand horizontal dashed lines identify the CEP shown in other scans (φ = π/2). Lines
are drawn for delayed resonances as outlined in Table 3.1.

. . . . . . . . . . . . . . 85

3.17 HHG spectra as a function of (top) peak intensity at λ = 800 nm and (bottom)
driving wavelength at I = 3.8 × 1013 W/cm2 , in atomic hydrogen for a 20 o.c. pulse.
We applied a Blackman window to dipole acceleration, and results for each intensity(top) or each wavelength (bottom) have been normalized to the maximum signal
in the region between 8.75ω and 9.5ω. Overlaid curves denote the ∆ω prediction,
Eq. (3.2), for (top) delayed resonances as outlined in Table 3.1 and (bottom) 9ω → 3p
resonances delayed by 0.5 (solid black), 1.8 (dashed black), and 3.0 (dotted black)
optical cycles.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

xviii
3.18 Temporal profile of radiation within 1ω of the 9th harmonic for a selection of electricfield envelopes (horizontal bars) in atomic hydrogen for an 800 nm, 20 o.c. pulse at
3.8 × 1013 W/cm2 . We note that the results have not been normalized. A verticaldashed line denotes the peak of the laser. The center bar corresponds to the sinesquared envelope used elsewhere in this thesis. In the lower panel, lines are drawn
for delayed resonances as outlined in Table 3.1.

. . . . . . . . . . . . . . . . . . . . 88

3.19 Temporal profile of radiation within 1ω of the 9th harmonic as a function of atomic
nuclear charge Z in a hydrogen-like atom while also varying wavelength (λ = 800Z2nm) and intensity (I = Z 6 ∗ 3.8 × 1013 W/cm2 ) to keep resonance times fixed. A
20 o.c. pulse was used throughout. For each Z the results have been normalized
to a maximum of 1 for sake of comparison (note: spectral intensity here scales√
roughly ∝ Z). A vertical dashed line denotes the peak of the laser and a horizontal
dashed line identifies the pulse parameters used as a base in this Chapter (i.e.,
Z = 1, I = 3.8 × 1013 W/cm2 , λ = 800 nm). Overlaid with lines denoting resonancetimes for 7- (magenta) and 8- or 9- (white) photon resonances. . . . . . . . . . . . . 92
3.20 Temporal profile of radiation within 1ω of the 9th harmonic for scans of: (top) atomic
nuclear charge Z and wavelength λ, (middle) Z and intensity I, and (bottom) λ andI. For all scans, the 9ω → 3p resonance time is kept fixed. A 20 o.c. pulse wasused throughout. For each element of the scan (i.e., horizontal slice), the results
have been normalized to a maximum of 1 for sake of comparison. A vertical dashedline denotes the peak of the laser and a horizontal dashed line identifies the pulse
parameters used as a base in this Chapter (i.e., Z = 1, I = 3.8×1013 W/cm2 , λ = 800
nm). Lines are drawn for delayed resonances as outlined in Table 3.1. . . . . . . . . 94
3.21 Same as top panel of Fig. 3.20, but with no normalization for each Z and on alogarithmic scale. Unscaled versions for other panels show similar trends, except
with an overall ≈ I 4 scaling of the radiation magnitude. . . . . . . . . . . . . . . . . 95

xix
3.22 Temporal profile of radiation within 1ω of the 9th harmonic for scans of pulse duration, Nτ , in optical cycles and (top) atomic nuclear charge Z, (middle) wavelength
λ, and (bottom) intensity I. For all scans, the 9ω → 3p resonance time is kept fixed.
For each element of the scan (i.e., horizontal slice), the results have been normalized
to a maximum of 1 for sake of comparison. A vertical dashed line denotes the peakof the laser and a horizontal dashed line identifies the pulse parameters used as abase in this Chapter (i.e., Z = 1, Nτ = 20 o.c., I = 3.8 × 1013 W/cm2 , λ = 800 nm).
Lines are drawn for delayed resonances as outlined in Table 3.1.

. . . . . . . . . . . 98

3.23 Same as middle panel of Fig. 3.22, but with no normalization for each Nτ and on
a logarithmic scale. Unscaled versions for other panels show similar trends, except
with an overall ≈ I 4 scaling of the radiation magnitude. . . . . . . . . . . . . . . . . 99
3.24 Temporal profile of radiation within 1ω of the 9th harmonic as a function of driving
pulse CEP in atomic hydrogen for an 800 nm, 4 o.c. pulse at 3.8 × 1013 W/cm2 .
A vertical dashed line denotes the peak of the laser. Lines are drawn for delayedresonances as outlined in Table 3.1.

. . . . . . . . . . . . . . . . . . . . . . . . . . . 99

3.25 Same as Fig. 3.5, but extending to durations up to 50 o.c..

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

3.26 High harmonic spectrum generated in helium, exhibiting a rich structure in the nearthreshold region of the spectrum between harmonics 5 and 12. The helium ionizationthreshold of 0.94 a.u. = 8.1ω, where ω = 0.114 a.u. corresponds to the central driverwavelength of 400 nm. Other laser parameters: sin2 pulse of 20 optical cycles withpeak intensity of 3.3 × 1014 W/cm2 . Note that the intensity of the spectrum is scaledby the signal at the fundamental frequency. (cf. hydrogen results in Fig. 3.1).

. . . 102

xx
3.27 HHG spectra as (top) a function of driving wavelength at I = 3.3 × 1014 W/cm2 and(bottom) a function of peak intensity at λ = 400 nm. Nτ = 20 o.c. for both panels.
Results for each wavelength (top) or each intensity (bottom) have been normalized
to the maximum signal in the region between 8.75ω and 9.5ω. Red dashed lines
denote field-free excited state transition energies. (cf. hydrogen results in Fig. 3.2
and Fig. 3.3).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

3.28 Temporal profile of radiation within 1ω of the 9th harmonic as a function of (top)pulse duration, (middle) peak intensity, and (bottom) driving wavelength. All inhelium atom for a pulse of (unless specified) 20 o.c., 400 nm, and 3.3 × 1014 W/cm2 .
For each element of the scans (i.e., horizontal slice), the result has been normalized
to a maximum of 1. A vertical dashed line denotes the peak of the laser and ahorizontal dashed line identifies the base pulse parameters.

. . . . . . . . . . . . . . 104

3.29 High harmonic spectra generated in atomic hydrogen. (top) Broad plateau extending
to cut-offs (near 27 and 67ω, denoted by dashed lines) with many peaks appearing
offset from N ω; (bottom) near-threshold region exhibiting a rich structure. The
hydrogen ionization threshold of 0.5 a.u. = 17.5ω, where ω = 0.028 a.u. corresponds
to the central driver wavelength of 1600 nm. Other laser parameters: sin2 pulse of
10 optical cycles with peak intensities of 0.75 × 1013 and 4.95 × 1013 W/cm2 . Notethat the intensity of the spectra are scaled by the signal at the fundamental frequency. 106
3.30 Temporal profile of radiation within 1ω of the (top) 13th and (bottom) 17th harmonics as a function of driving pulse intensity in atomic hydrogen for a 1600 nm,
10 o.c. pulse. For each intensity the results have been normalized to a maximum of
1 for sake of comparison. A vertical dashed line denotes the peak of the laser andhorizontal dashed lines identify the intensities in Fig. 3.29.

. . . . . . . . . . . . . . 107

xxi
3.31 Temporal profile of radiation within 1ω of the 13th harmonic as a function of driving
pulse intensity in atomic hydrogen for a 1600 nm, 10 o.c. pulse. For each intensity the
results have been normalized to a maximum of 1 for sake of comparison. A verticaldashed line denotes the peak of the laser and a horizontal dashed line identifies thelower intensity in Fig. 3.29 (i.e., 0.75×1013 W/cm2 ). Panels show identical radiation
data, but the lower is overlaid with lines denoting resonance times for 12−14ω (white)and > 14ω (magenta) resonances (there are no excluded resonances with a smalleror larger number of photons involved for these parameters) for 2s − 6p.
4.1. . . . . . . 108

Example of single-atom high harmonic spectrum generated in atomic hydrogen, obtained via numerical solution of the time-dependent Schrödinger equation, exhibitinga rich structure in the near-threshold region of the spectrum between harmonics 5
and 12. The hydrogen ionization threshold of 0.5 a.u. = 8.77ω, where ω = 0.057 a.u.
corresponds to the central driver wavelength of 800 nm. Other laser parameters:
sin2 pulse of 20 optical cycles with peak intensity of 7 × 1013 W/cm2 . Note that theintensity of the spectrum (as well as other spectra in this Chapter) is scaled by thesignal at the fundamental frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

4.2(a) Comparisons of results of interpolation (dashed lines) and exact numerical results(solid lines) for high harmonic spectrum and (b) for dipole acceleration, and (c) error
in harmonic phase, generated by a sin2 pulse of 20 optical cycles at wavelength of
800 nm and peak intensity of 9.1431 × 1013 W/cm2 . Interpolation performed on datawith intensities of 0.1, 0.2, 0.3, ..., 10.0 × 1013 W/cm2 . ((a) taken from [6].). . . . . 119

xxii
4.3

Trend of maximum absolute error in interpolated high harmonic spectrum (as compared to numerical TDSE result) as a function of the intensity sampling spacing dI(i.e., sampled intensities dI, 2dI, . . . , ∼ 1014 W/cm2 ). Shown are exemplary results
at peak intensities 3.97 × 1013 W/cm2 (blue circles with dashed line) and 9.14 × 1013
W/cm2 (red circles with dotted line). Other laser parameters as in Fig. 4.2. Notethat at sufficiently small intensity spacings, the error roughly drops as O(dI 4 ) (black
dashed lines).

4.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

Spectra for a selection of sample CEPs, generated by a sin2 pulse of 20 optical cycles
at wavelength of 800 nm and peak intensity of (top) 3.8 × 1013 W/cm2 and (bottom)
7.0 × 1013 W/cm2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

4.5

Spectral phase difference for a selection of sample CEPs, generated by a sin2 pulseof 20 optical cycles at wavelength of 800 nm and peak intensity of (top) 3.8 × 1013
W/cm2 and (bottom) 7.0 × 1013 W/cm2 . (Top panel taken from [6].)4.6. . . . . . . . 124(top) Spectrum and (bottom) spectral phase difference for a selection of sample
CEPs, generated by a short sin2 pulse of 4 optical cycles at wavelength of 800 nm
and peak intensity of 3.8 × 1013 W/cm2 .

4.7. . . . . . . . . . . . . . . . . . . . . . . . 126

Angle-resolved macroscopic high harmonic spectra from uniformly-separated hydrogen atoms in different configurations: (a) row of 100 atoms along the propagation
direction, (b) row of 100 atoms perpendicular to the propagation direction, and (c)lattice grid of 100 × 100 atoms in the plane spanned by the combination of (a) and(b). A Gaussian beam (beam waist 30 µm) with 20-cycle duration at 800 nm and9 × 1013 W/cm2 and an atom separation of L = 300 nm were considered. The red
and white lines represent the constructive interference conditions for atoms locatedalong and perpendicular to the propagation direction, respectively (see Eqs. (4.11)
and (4.10)). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

4.8

Comparison of spectra taken from the results in Fig. 4.7(c) along θd = 0.65 rad (bluesolid line) and in forward direction (θd = 0 rad, black dashed line).

. . . . . . . . . 130

xxiii
4.9

Angle-resolved macroscopic radiation spectrum for a gas jet of hydrogen atoms driven
by an 800 nm laser with peak intensity 7.0 × 1013 W/cm2 for different numbers ofcells (rows) and cell sizes (columns): (top) 103 cells, (middle) 104 cells, and (bottom)
105 cells. (left) point-like radiators, (middle) 3 µm cells, and (right) 5 µm cells.

. . 134

4.10 Convergence of macroscopic on-axis radiation with number of cells (cell width has
no effect), generated by a sin2 pulse of 20 optical cycles at wavelength of 800 nm andpeak intensity of (top) 3.8 × 1013 W/cm2 and (bottom) 7.0 × 1013 W/cm2 . Verticaldashed lines denote the ionization potential (left line) and HHG cut-off (right line).

135

4.11 Comparison of single-atom TDSE results and macroscopic TDSE interpolation results for high harmonic generation in hydrogen atoms. On the left, the (a) singleatom and (b) macroscopic on-axis spectrum is shown as a function of peak intensity.
The spectra on the right show the comparison (single atom: black dashed line,
macroscopic: red solid line) at peak intensities (c) I0 = 3.8 × 1013 W/cm2 and (d)I0 = 7.0 × 1013 W/cm2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
4.12 Comparison of macroscopic calculations in which the single-atom TDSE results areadded coherently (solid red line) and incoherently (dot-dashed blue line) for highorder harmonic generation at I0 = 3.8 × 1013 W/cm2 . . . . . . . . . . . . . . . . . . 138
4.13 Angle-resolved near-threshold part of the macroscopic radiation spectrum for a gas
jet of hydrogen atoms driven by an 800 nm laser with peak intensity 3.8 × 1013
W/cm2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.14 Wavelet analysis of the macroscopic radiation for a gas jet of hydrogen atoms driven
by an 800 nm laser with peak intensity 7 × 1013 W/cm2 at detection angles (a)θ = 0 mrad (on-axis), and (b) θ = 5 mrad. The shifted energy of the 2p, 3p, and4p states relative to the 1s (ground) state are denoted by red lines, while the region
containing the higher Rydberg states (shifted 5p to threshold) are shaded in red.
Wavelet parameter Ω0 = 54.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

xxiv
4.15 Comparison of (top) single-atom results and (bottom) on-axis macroscopic TDSE
interpolation results for near-threshold HHG spectra as a function of peak intensity
for an 800 nm, 20 o.c. pulse interacting with atomic hydrogen. Results for eachintensity have been normalized to the maximum signal in the region between 8.75ω
and 9.5ω. Horizontal red dashed lines denote field-free excited state transition energies and vertical dashed lines identify intensities selected in Sec. 3.4: 2.0 × 1013 ,3.0 × 1013 , 3.8 × 1013 , 5.6 × 1013 , and 7.0 × 1013 W/cm2 .

. . . . . . . . . . . . . . . 142

4.16 Comparison of (top) single-atom results and (bottom) on-axis macroscopic TDSE
interpolation results for temporal profile of radiation within 1ω of the 9th harmonic
as a function of driving pulse intensity in hydrogen for an 800 nm, 20 o.c. pulse.
For each intensity, the results have been normalized to a maximum of 1 for sakeof comparison. A vertical dashed line denotes the peak of the laser and horizontal
dashed lines identify intensities selected in Sec. 3.4: 2.0 × 1013 , 3.0 × 1013 , 3.8 × 1013 ,5.6×1013 , and 7.0×1013 W/cm2 . Lines are drawn for delayed resonances as outlinedin Table 3.1.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143

4.17 Spatial distribution of carrier-envelope phase in a focused broadband Gaussian laser
pulse: (a) Gouy phase (g0 = 0) and (b) focal phase with g0 = −2. For the latter, reference lines are drawn showing (magenta) position where laser intensity falls below
1% of peak, (cyan) position where gas density ρ falls below 10% of peak for a gas
width of 0.23z0 , and (white) estimates of the effective slit locations based on interference patterns (0.925 times the transverse node separation near z = 0 in Eq. (4.12)).
(Taken from [6].). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

4.18 Comparison of simulated angle-resolved macroscopic high harmonic spectra using(a) Gouy phase (g0 = 0) and (b) focal phase with g0 = −2. Lines denote peaks√of a double-slit pattern with slit separation of 0.925 2W0 ≈ 39µm in the direction
perpendicular to the laser propagation. Laser and gas jet parameters are given in
the text. (Taken from [6].). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147

xxv
4.19 Position of first minimum (circles) and first maximum (squares) in the interferenceof the 7th harmonic as a function of the (a) width of the central gas jet σz , (b)g0 , and (c) W0 . Other parameters are kept the same as in Fig. 4.18. The line in
the bottom panel shows the prediction for the first side maximum based on the
double-slit formula with estimated slit separation of 0.925L0 . (Taken from [6].) . . . 149
4.20 Comparison of single-atom results (black dotted) and macroscopic TDSE interpolation results for different driving laser CEP profiles: constant phase (purple), Gouy
phase (red), focal phase with g0 = −2 (green). For a gas jet of hydrogen atoms
driven by an 800 nm laser with peak intensity 7.0 × 1013 W/cm2 .

. . . . . . . . . . 151

4.21 Convergence of macroscopic on-axis radiation with number of cells for (top) Gouy
phase and (bottom) focal phase with g0 = −2, generated by a sin2 pulse of 20 optical
cycles at wavelength of 800 nm and peak intensity of 7.0 × 1013 W/cm2 . Verticaldashed lines denote the ionization potential (left line) and HHG cut-off (right line).
Note different cell numbers than Fig. 4.10.

. . . . . . . . . . . . . . . . . . . . . . . 152

4.22 Comparison of (dot-dash black) single-atom results versus (dashed blue) incoherentand (solid red) coherent macroscopic summation, for HHG in a gas jet of hydrogenatoms driven by a 1600 nm laser pulse of 10 optical cycles with peak intensity
8.8 × 1013 W/cm2 . Note that the intensity of the spectra are scaled by the signal at
the fundamental frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
4.23 Macroscopic radiation spectra for a gas jet of hydrogen atoms driven by a 1600
nm laser. (top) Angle-resolved spectrum for a peak intensity of 8.8 × 1013 W/cm2 .
(bottom) On-axis spectra shown as a function of peak intensity.

. . . . . . . . . . . 156

xxvi
4.24 Comparison of single-atom TDSE results and macroscopic TDSE interpolation results for high harmonic generation in helium atoms. On the left, the (a) singleatom and (b) macroscopic on-axis spectrum is shown as a function of peak intensity.
The spectra on the right show the comparison (single atom: black dashed line,
macroscopic: red solid line) at peak intensities (c) I0 = 3.3 × 1014 W/cm2 and (d)I0 = 8.0 × 1013 W/cm2 . (cf. hydrogen results in Fig. 4.11). . . . . . . . . . . . . . 157

4.25 Angle-resolved near-threshold part of the macroscopic radiation spectrum for a gas
jet of helium atoms driven by a 400 nm laser with peak intensity 3.3 × 1014 W/cm2 .
(cf. hydrogen results in Fig. 4.13). . . . . . . . . . . . . . . . . . . . . . . . . . . . 158

4.26 Comparison of (top) single-atom results and (bottom) on-axis macroscopic TDSE
interpolation results for near-threshold HHG spectra as a function of peak intensity
for a 400 nm, 20 o.c. pulse interacting with helium. Results for each intensity
have been normalized to the maximum signal in the region between 8.75ω and 9.5ω.
Red dashed lines denote field-free excited state transition energies, and vertical lines
identify select intensities: 2.4×1014 , 3.3×1014 , and 8.0×1014 W/cm2 . (cf. hydrogen
results in Fig. 4.15). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160

4.27 Comparison of (top) single-atom results and (bottom) on-axis macroscopic TDSE
interpolation results for temporal profile of radiation within 1ω of the 9th harmonic
as a function of driving pulse intensity in helium for a 400 nm, 20 o.c. pulse. For
each intensity, the results have been normalized to a maximum of 1 for sake ofcomparison. A vertical dashed line denotes the peak of the laser and horizontal
dashed lines identify select intensities: 2.4 × 1014 , 3.3 × 1014 , and 8.0 × 1014 W/cm2 .
(cf. hydrogen results in Fig. 4.16). . . . . . . . . . . . . . . . . . . . . . . . . . . . 161

Chapter 1
Introduction

1.1

Background
Lasers are widely used in medicine, industry, defense, and scientific research. The ability to

precisely apply power is important itself, but lasers offer other benefits, especially in the broad realmof detection and measurement. Observing reflection of light off of and absorption into materials has
consistently exposed otherwise inaccessible information. One form of light is particularly relevant
for its uses in this way: x-ray imaging and x-ray diffraction off crystals. The short wavelength andhigh photon energy of x-ray radiation makes it a tool for detailed observation and analysis thateasily pierces obscuring materials (e.g., flesh). Ultraviolet (UV) lasers tuned to energies relevantin atoms and molecules serve as a different sort of probe. The main interest of this thesis is inmodeling the generation of ultrashort UV laser pulses from a theoretical perspective with the goalof deepening our understanding and control of the sources and the underlying physical processes.

1.1.1

Single- and multi-photon absorption

In 1887, it was observed that metals exposed to light would release electrons with kinetic
energy related only to the color (frequency) of the light and not its intensity [7]. In 1905, Albert
Einstein described this effect as resulting from the quantization of light [8]. Electrons would absorb
energy from light in the form of single photons with energy proportional to the frequency, butwould not absorb more than one photon no matter the photon density (i.e., intensity). In 1930,
Maria Göppert-Mayer suggested that a sufficiently intense light would allow electrons to absorb

2
multiple photons, reaching higher kinetic energies than allowed by Einstein’s rule [9]. With theadvent of lasers in 1961, her prediction was verified – two-photon absorption was observed [10]. As
laser technology advanced to allow higher intensities, multiphoton absorption involving more than
two photons became possible.

1.1.2

Applications

Short duration, high energy laser pulses have applications in many disciplines [11, 12]. First,
lasers with shorter wavelength (higher frequency) allow use of light as “microscopes” with finer
spatial resolution. Likewise, taking “pictures” of rapid processes requires the use of measurement
techniques that operate on a similar time scale, or the behavior of interest will be blurred out.
Finally, controlling a process requires tools that can be applied at the characteristic frequency orenergy scale of the process. X-rays in general have wavelengths smaller than 10 nm and energies
above 100 eV; the subject of this thesis, high-order harmonic generation, in particular leads topulses with durations in the attosecond (10−18 s) to femtosecond (10−15 s) range [13, 14]. This
time scale corresponds to that of single-electron motion in valence shells and bulk-electron motion inlarger molecules – such as those found in biology. In combination with fine spatial resolution, HHG
lasers can benefit chemical and biological research and engineering via their ability to “film” andcontrol chemical reactions [15, 16]. Other methods of creating x-ray lasers, involving kilometerlong electron accelerators, have greatly increased our understanding of proteins, photosynthesis,
and many other chemical reactions. In contrast, HHG can produce well-controlled lasers – witha wide range of wavelengths – on a tabletop scale, allowing scientists in smaller labs access to a
powerful research tool. Finally, the process of HHG itself can be used to probe atomic and molecular
structure [17–19].

1.1.3

Hartree atomic units

When discussing physical phenomena, it is useful to use a unit system with an appropriate
scale. For example, consider various temperature systems – Fahrenheit is useful for gauging humanDimension
ChargeMassEnergyAngular momentum
LengthTime
Permittivity
Speed of light

Constant
electron charge, eelectron mass, me
2x hydrogen ionization potential~
Bohr radius, a0
1st Bohr orbital period / 2π
Coulomb’s constant ke = 1/4πε0c = 1/α = 137.03500

3

Conversion (1 a.u. = ... )1.60218 × 10−19 C
9.10938 × 10−31 kg
27.2110 eV
1.05457 × 10−34 J*s
0.05292 nm
24.18884 as
8.98755 × 109 Nm2 /C2Table 1.1: Description of Hartree atomic units, including reference constants and conversion tocommon units.

comfort and safety (corresponding roughly to 0◦ − 100◦ F ), while Celsius more intuitively describes
water (freezing and boiling at 0◦ C and 100◦ C, respectively), and Kelvin is important for applications
where total thermal energy is relevant (since temperatures in K are proportional to the thermalenergy, unlike in Fahrenheit or Celsius).
In strong-field atomic physics, one natural scale is that of the electron. Rather than talking
about masses smaller than 10−25 kilograms, charges as multiples of 1.6 × 10−19 C, and so forth, we
use a unit system where as many dimensions as possible are scaled so that electronic behavior isdescribed in numbers of order 1. Most relevant constants are set to 1, but consistency is maintained,
requiring some (notably the speed of light, c) to have other values. Generally, constants with unity
values are left out of equations. Also, the unit abbreviation, a.u., is the same for all atomic units
regardless of dimension. In certain circumstances, equivalence is sensible (e.g., a photon withfrequency ω a.u. will always also have energy ω a.u.), but in other cases is not (that same photon
will have wavelength λ = 2πc
ω a.u.). Table 1.1 includes a few of the most important conversions.

1.2

Strong field ultrafast physics
Classically, a free electron in a laser field of the form
E(t) = E0 sin(ωt),

(1.1)

4
with intensity I = E02 and frequency ω, will move as follows
qE= −E0 sin(ωt)
Zm
E0
cos(ωt) + v0
v =
adt =
ω
Z
E0
x =
vdt = 2 sin(ωt) + v0 t + x0
ω
a =

(1.2)
(1.3)
(1.4)

where x0 and v0 are the initial position and velocity, respectively. From this, we can define characteristic values for the field. First, the characteristic length of electron motion is simply the
amplitude of the oscillation – the so-called quiver radius
rq =

E0
.
ω2

(1.5)

Similarly, we can consider a characteristic kinetic energy associated with field motion, the ponderomotive energy UP . Excluding energy associated with the average velocity of the electron (or,transforming to a frame where this is zero), we have
T
UP

1
E2
mv 2 = 02 cos2 (ωt),
2
2ω
E02
I
= hT i =
=
.
4ω 2
4ω 2
=

(1.6)
(1.7)

The ponderomotive energy in particular appears often, as well as being used to define several
dimensionless parameters [20]:
r
γ =
zf

=

z =

IP
2UP
UP
2mc2
UP
,
ω

(1.8)
(1.9)
(1.10)

where the Keldysh parameter γ describes ionization as occurring through tunneling (γ  1) ormultiphoton ionization (γ  1); the free-electron intensity parameter zf describes the significanceof relativistic effects (minimal for all fields used in this thesis, where zf  10−3 ); and, the nonperturbative intensity parameter z naturally describes how non-perturbative interactions will be.
z will also appear later in its role as a measure of the typical number of photons involved in theinteraction with the field.

5
1.2.1

Laser pulses

In contrast to the idealized continuous wave (CW) lasers above (Eq. (1.1)), strong-field lasers
have a finite extent in time and driving laser pulses are given by:
E(t) = f (t)E0 sin(ωt + φ),
where E0 =√(1.11)

I is the amplitude of the electric field for peak intensity I, ω = 2πc/λ is the central

frequency for wavelength λ, φ is the carrier-envelope phase (CEP), and f (t) is some envelope
function, such as the sine-squared or Gaussian envelope2

fsin2 (t) = sin



πt
τ


(1.12)
fGauss (t) = exp − ln(2)

2t
τ

2 !
(1.13)

with τ = 2πNτ /ω being the pulse duration for a pulse of Nτ cycles. As is typical, we define
durations as full width for sine-squared, and full width half max (FWHM) for Gaussian envelopes.
In this thesis, we do not report results with Gaussian envelopes, but do additionally use variationsof the sine envelope with faster rise or fall:

fsinn sinm (t) =





sinn

πt
τ




sinm

πt
τ



t≤0
(1.14)

t≥0for even integers n, m.

1.2.2

Tunnel and multiphoton ionization

With laser pulses of moderate intensity (roughly 1013 −1015 W/cm2 ) and wavelengths near thevisible (hundreds of nanometers), two pathways contribute to ionization, with dominance predictedroughly by the Keldysh parameter γ: multiphoton ionization (γ > 1, at shorter λ and lower I) andtunnel ionization (γ > 1, at longer λ and higher I). The diagram in Fig. 1.1 shows both processes.
An electron is initially in the ground (1s) state of an atomic potential (black curve).

6
In the case of tunnel ionization, the potential is distorted (blue curve), leading to a finite
potential barrier. After tunnel ionization, the electron is assumed to start at rest at the tunnel exit,
x0 . The rate of ionization depends exponentially on the instantaneous value of the electric field[21], resulting in bursts of ionization every half-cycle. In the case of multiphoton ionization (threephoton ionization shown here), the electron absorbs several photons of energy ω (purple arrows),
resulting in non-zero kinetic energy, T0 at the moment of ionization. For weak fields, perturbation
theory gives a rate of I Nγ , where Nγ is the number of photons involved in ionization. However,
non-perturbative, above-threshold ionization (ATI) – in which more photons than the minimum
required for ionization are absorbed – is an area of active research not discussed in this thesis (cf.
[22–24]).

1.3

AC Stark shift1
In the presence of an alternating electric field, the energies of states shift. For states whose

eigenenergy En is much less than the photon energy ω, this shift becomes simply UP [25]En0 = En + UP .

(1.15)As a result, the effective ionization potential also shifts
IP0 = IP + UP ,(1.16)

and as a consequence, multiphoton ionization channels may close. In Fig. 1.2, we show the shifted
energies (in ω) of ionization (IP0 , black dashed line) and excitation from 1s to hydrogen n = 2
0 − E and E 0 − E , blue line) and n = 3 (e.g., E 0 − E , orange line) shells. When the(i.e., E2s
0
0
0
2p
3s

minimum energy to ionize an electron exceeds the energy of an Nγ -photon process (i.e., IP0 > Nγ ω),multiphoton ionization must occur primarily through absorption of Nγ + 1 (or more) photons, andwe say the Nγ -photon channel has closed. For transition into bound states, resonant excitation
has greatly enhanced probability if En0 − E0 = Nγ ω. For short intense pulses, the energy shift as
1

The results of this Section are also presented in [3].

7Figure 1.1: Diagram of ionization processes. (black) Field-free atomic potential and 1s bound stateenergy. (blue) Tunnel ionization: Potential distorted by external field, resulting in tunnel exit x0 .
(purple) Three-photon ionization resulting in excess kinetic energy T0 .

8a function of time may also be relevant (Fig. 1.3). For this purpose, we define an approximate
time-dependent shift [26, 27]
En0 (t) = En + UP (t)(1.17)

where
UP (t) ≡2~|E(t)|.
4ω 2(1.18)

In [3], we investigated the impact of such channel closings and resonances on the distributionof Rydberg state population after a pulse. To this end, we considered peak intensities such thatn = 8 states shifted into multiphoton resonance with Nγ = 10, 11, . . . 15 at central frequency
ω, corresponding to a wavelength of 800 nm. Due to the finite duration of laser pulses (versus
monochromatic CW fields), more than one manifold of states can be resonantly excited. For
example, the bandwidth of a 20 o.c. sin2 pulse at 800 nm covers all excited states n ' 6 within thesame Nγ photon process. We further note that in the present cases, the n = 3 and n = 2 states areapproximately resonant via Nγ − 1 and Nγ − 2 photon processes, respectively, assuming that theenergy shift of these states equals the ponderomotive energy as well.
The population in each quantum state (n, l) at the end of the pulse is shown in Fig. 1.4 forintensities between 1013 and 1014 W/cm2 . The results shown in the left (right) column correspondto cases in which the Rydberg states are resonant with an even (odd) number of photons. In
all cases, we see that the highest angular momentum states (l > 7) are not much populated, inagreement with previous studies [28] and semi-classical estimates [29]. In the results, we observe
signatures of selection rules resulting in dominant population of states with an even (odd) angularmomentum quantum number for the absorption of even (odd) photons in the plots on the left
(right). This shows that in the intensity regime of 1013 to 1014 W/cm2 the parity selection rules,
previously studied for monochromatic fields and flat-top pulses, are effective for long pulses of about
20 cycles. However, we observe that the contrast between the population in even and odd states isstronger at lower intensities.
The results also show that for an odd parity process (right column) predominantly one angular9

Figure 1.2: Diagram of the AC Stark shift as a function of laser intensity for a wavelength of 800
nm. (solid lines) Excitation and (dashed line) ionization energies are given in units of ω. As theintensity increases, N -photon ionization channels close (e.g., the 9-photon channel at intensities
above about 6 × 1012 W/cm2 ) and excited states shift across resonances.

Figure 1.3: Diagram of the time-dependent AC Stark shift. Near the peak of the pulse, the effective
ionization potential (top dashed line) rises above the energy of five photons (red arrows), while the2s state energy (top solid line) shifts briefly into 4-photon resonance.

10
momentum state (l = 5) is populated. This is in agreement with the results presented by Li etal. [30], who conjectured that electrons in the low angular momentum states more easily absorb
additional photons resulting in suppression of population in these states due to ionization. However,
our results for even parity processes exhibit a pattern, alternating in l, showing that both low andhigh angular momentum states, except for the s-states, remain populated at the end of the pulse.
Chapter 3 investigates related processes in resonance-enhanced HHG.

1.3.1

High-order harmonic generation

In high-order harmonic generation (HHG), an atom or molecule is exposed to a high-intensity,
low-frequency laser. An electron absorbs many photons (up to thousands [31]) before being driven
back into the atom; all of the absorbed energy is then emitted as a single photon. When a laser
shines on a gas consisting of many such atoms, the emitted photons form high-frequency laser lightwith lower intensity. The emitted photons have a large range of frequencies with similar intensities– an unusual and useful trait. Of primary interest are the high-energy portions (into the x-ray
regime) and the incredibly short duration of these laser pulses.
HHG is characterized by a so-called plateau where radiation occurs with consistent intensity at odd multiples of the driving laser frequency, ω. In Fig. 1.5, we show a schematic HHG
spectrum. The first few harmonics (1ω, 3ω, 5ω, and 7ω) below the ionization threshold (IP ) show
exponentially-decreasing intensity as predicted by perturbation theory. However, the next several
odd harmonics (through 21ω) appear at the same intensity, forming the plateau. After a so-called“cut-off” energy (21ω), the intensity of harmonics again decreases exponentially.
HHG can be understood through a three-step semi-classical model [32–34] (Fig. 1.6). An
electron is initially in the ground state of an atomic Coulomb potential. Then, an intense fieldinitiates tunnel ionization. The freed electron oscillates purely under the influence of the external
field (Eq. (1.4)), gaining kinetic energy. When the electron again approaches the vicinity of thenucleus, it has a chance to recombine back to the ground state. In doing so, it emits the excessenergy as a single high-energy photon. Semi-classical trajectory analysis (Sec. 1.4.1) predicts the

11

14
13
12
11
10
9
8
7
6
5
4
3
2
1

10−5

n

n

10−6

14
13
12
11
10
9
8
7
6
5
4
3
2
1

10−7

l

10−5

0
1
2
3
4
5
6
7
8
9
10
11
12
13

(c)

l

10−6

l

10−3

10−4

10−5

(d)
0
1
2
3
4
5
6
7
8
9
10
11
12
13

n

10−4

14
13
12
11
10
9
8
7
6
5
4
3
2
1

10−7

n

14
13
12
11
10
9
8
7
6
5
4
3
2
1

(b)
0
1
2
3
4
5
6
7
8
9
10
11
12
13

10−8

0
1
2
3
4
5
6
7
8
9
10
11
12
13

(a)

10−6

l

Figure 1.4: Excited state distribution as function of quantum numbers n (vertical axis) andl (horizontal axis) at the end of 800 nm, 20 cycle pulses with sine-squared envelope and peakintensities: (a) I0 = 3.4 × 1013 W/cm2 , (b) I0 = 6.0 × 1013 W/cm2 , (c) I0 = 8.6 × 1013 W/cm2 ,and (d) I0 = 11.2 × 1013 W/cm2 . Left (right) column corresponds to cases in which the Rydberg
states are resonant with an even (odd) number of photons. (Taken from [3].)12Figure 1.5: Diagram of a high-order harmonic generation spectrum marked by a plateau of radiation
at similar intensities.

13

Figure 1.6: Diagram of the 3-step model of high-order harmonic generation. (Taken from [4].)maximum kinetic energy upon return as 3.17UP , resulting in the radiation cut-off of IP + 3.17UP .
Various additional considerations result in a more complicated interaction and may allow for higher
return energies, including inclusion of the Coulomb force in classical trajectories, the presence ofquantum trajectories, and multiphoton ionization.
As with perturbative harmonics, atomic inversion symmetry across the laser polarization
direction prohibits even-harmonic radiation [35]. This symmetry is not present in other systems(e.g., most molecules) or when significant electron density departs the ground state.
HHG has been a subject of intense research for decades, for a variety of reasons. Firstly, unlike
other methods of generating high energy photons (e.g., linear accelerators such as at SLAC [36]),
HHG can be performed on the tabletop scale. Additionally, HHG allows a higher level of coherence
and control of various parameters including the CEP, polarization, and angular momentum of thelaser [37, 38].

14
1.3.2

Free-induction decay

In the context of strong-field physics, free-induction decay (FID) is simply the process
whereby electrons in excited np states decay to the ground state and emit photons with the corresponding energy ω = En − E0 [39]. The bandwidth of this radiation is generally quite small relative
to that of harmonic radiation, appearing as sharp peaks. Additionally, FID continues well after
the end of the driving laser pulse, unlike most other radiative processes. In our numeric methods,
FID is not fully captured since the vacuum field necessary for spontaneous emission is absent –FID radiation appears due to the oscillation of electrons in a superposition state, but the excited
population never decays to the ground state and energy is not conserved. Fortunately, the timescaleof spontaneous decay is generally orders of magnitude longer than the duration of our pulses, andthe approximation of infinite lifetime is good.

1.4

Numeric methods
This thesis does not cover the development of new methods for simulation of single-electron

behavior, but we briefly describe tools developed by others – in some cases implemented during thecourse of this research.

1.4.1

Semi-classical trajectories

From the 3-step model of HHG, one may derive an approximative description where electron
motion is simulated through simple Newtonian methods. Namely, the initial position and momentum of the electron are set by ionization conditions (typically, placing the electron at the tunnel
exit with zero initial velocity) and motion then proceeds according to the classical equations of motion, often excluding the effect of the Coulomb potential. In this case, motion may be analytically
predicted (i.e., Eq. (1.4)), but when the full force is included:
F (x, t) = FE (t) + FC (x) = −E0 sin(ωt) −

x
,
|x3 |

(1.19)

15
numeric solutions are required, such as Verlet integration [40] where the position and velocity arecomputed for progressive steps in time, ti+1 = ti + dt:
x(ti+1 ) = x(ti ) + v(ti )dt + 0.5a(ti )dt2

(1.20)

a(ti+1 ) = F (x(ti+1 ), ti+1 )(1.21)

v(ti+1 ) = v(ti ) + 0.5 (a(ti+1 ) + a(ti )) dt.

(1.22)

Radiation is then predicted whenever x crosses 0, and the resulting energy is given by the kinetic
energy T = 1/2v 2 plus IP .
A broader discussion of semi-classical trajectories can be found in Hostetter et al. [41],including: methods to avoid issues with the singularity in the Coulomb field and derive the radiation
phase computed from classical action; and, an alternate set of initial conditions meant to represent
multiphoton ionization. We do not present results of calculations of semi-classical trajectories inthis thesis, but any remarks are based on a reproduction of the algorithm as presented in [41].

1.4.2

Finite difference and eigenstate problems

For the vast majority of Schrödinger equation cases, there is no analytic solution to eventhe time-independent problem. Thus, we proceed with numerical solution through discretizationof space and time, for example the wavefunction ψ(x, t) → ψ(xi , tj ) with indices i, j = 0, 1, ....
In this thesis, all discretization occurs on uniformly-spaced grids and with uniform time steps, so
xi+1 = xi + ∆x and tj+1 = tj + ∆t.
Spatial derivatives are then computed through finite differences, for example the 2nd-orderapproximations:∂ψ(xi ) ≈∂x∂2
ψ(xi ) ≈∂x2−ψ(xi−1 ) + ψ(xi+1 )∆xψ(xi−1 ) − 2ψ(xi ) + ψ(xi+1 ),∆x2(1.23)(1.24)

16
which may be expressed as matrix equations, such as 2 −1  ψ(x0 ) 






−1 2 −1  ψ(x1 ) 





∂ 2 ψ(xi )
1 
..
.
.
.



.
.
.
.
=
.
.
.
.


∂x2
∆x2 









−1 2 −1  ψ(xN −1 )





−1 2
ψ(xN )(1.25)The time-independent Schrödinger equation (in one dimension)∂2
En ψn (x) = − 2 + V (x) ψ(x)∂x(1.26)

becomes an eigenvalue/eigenvector problem. Our early eigenstate solutions were based on a trickof the time-dependent numeric propagator (called Imaginary Time Propagation, see Sec. 1.4.3.2).
However, more recently we have utilized [42] the Arnoldi method of eigenvalue solutions [43].
While the TISE has spherical symmetry for atomic potentials, we often require eigenstates
as the initial conditions for further time-dependent interactions with laser fields without such symmetry. For all results of this thesis, linearly-polarized lasers are used and the electron density
remains in the subspace where magnetic/azimuthal quantum number m = 0. Thus, we perform
computations with discretization in cylindrical (ρi , zi ) or spherical (ri , Yl0 ) coordinates.
1.4.3

Time-dependent Schrödinger equation

Solutions to the TDSE (length gauge)


∂
∇2
~
i ψ(~r, t) = Hψ(~r, t) = −+ V (~r) + E(t) ∗ ~r ψ(~r, t)∂t
2(1.27)

naturally require propagation in time. Considering the propagator formulation
ψ(~r, t) = e−iH∆t ψ(~r, t),(1.28)

we may achieve numerical tractability through the Cayley form [44]e−iH∆t ≈1 − i ∆t
2 H
1 + i ∆t
2 H,

(1.29)

17
and the Crank-Nicolson method [45]:




∆t
∆t
1 + i H ψ(t + ∆t) = 1 − i H ψ(t).
22(1.30)The resulting problem in multiple dimensions may be performed either through approximating theeffect of the Hamiltonian along each dimension (e.g., H = Hx1 + Hx2 ):∆t∆t

e−iH∆t ≈ e−iHx1 2 e−iHx2 ∆t e−iHx1 2 ,(1.31)or through direct solution relying on modern sparse-matrix methods, such as the Generalized
Minimal Residual method [46]. The former is simpler to code, but the latter benefits from better
parallelization on large computational systems. Results in this thesis comprise both methods, which
we have previously shown to give consistent results [3].

1.4.3.1

Exterior complex scaling

One issue with simulations on discrete grids is that the grid cannot have infinite extent – weare in practice placing our interaction in a box. For weak fields, it may be possible to use numericalgrids large enough that no appreciable electron density reaches the edge, but this is not feasible forsimulations in the strong-field regime of interest. Instead, one must ensure that behavior of interest(e.g., HHG trajectories) fits within the extend of the grid, and appropriately handle departureof ionized electron density. The former is generally a matter of checking numeric convergence by0toscanning the box size, but we have found a radius of three times the quiver radius rq = E
ω2generally suffice for HHG simulations. For the escaping density, we use Exterior Complex Scaling(ECS) [47, 48], where the edges of the grid are rotated into complex space by an angle η = π/4:




ρ
ρ ≤ ρ0ρ (or r) =(1.32)ρ0 + (ρ − ρ0 )eiη ρ ≥ ρ0

z=





−z0 + (z + z0 )eiη





z







z0 + (z − z0 )eiη18z ≤ −z0−z0 ≤ z ≤ z0

(1.33)z ≥ z0 .

We apply this transformation only to the kinetic energy term of the Hamiltonian to avoid undesirable density-explosion from the oscillating electric field [49].

1.4.3.2

Imaginary Time Propagation

Given a numeric library for solution of the TDSE, Imaginary Time Propagation (ITP) is a
simple and popular method of generating initial states similar to the power method. An initial
wavefunction guess ψ(~r, t = 0) is made based on analytic solutions or random generation. The
wavefunction is then propagated in steps of imaginary time ψi (~r, t = −iN ∆t) and normalizing|ψ|2 = 1 until reaching numeric convergence. This functions by replacing eigenstate quantum
oscillation:
ψ(~r, t) =X

cn ψn (~r, t = 0)e−iEn t(1.34)

nwith exponential enhancement of the ground state (i.e., state with most negative energy):
ψ(~r, t = −iN ∆t) =X

cn ψn (~r, t = 0)e−En N ∆t .

(1.35)

n

The first excited state may be acquired by projecting out the ground state after each iteration (viaGram-Schmidt orthogonalization), and so forth for additional states.
ITP has the advantage of building on tools already developed for (real) time propagation,and generally converges rapidly for the first few states. However, it becomes inefficient for higher
excited states (when repeated orthogonalization is required) and may cease functioning altogether
due to limitations of machine precision. For this reason, methods better-grounded in linear algebra(e.g., the Arnoldi method) are preferable.

19
1.4.3.3

Simulating spectra

To derive the radiation spectrum from simulation of electron motion, we evaluate the dipole
acceleration a(t) using the Ehrenfest theorema(t) = h∇V i ,and compute the spectrum as its Fourier transform:
Z T
2
P (ω) =
a(t)e−iωt dt ,(1.36)(1.37)

0where T is the final time of the simulation. Often, we apply a window function as described inSec. 1.5.1.
The dipole acceleration can generally be extracted from the integration of the wavefunction
at each time step during simulation without dramatically reducing efficiency. If necessary, it may
be output at a lower time resolution (dt = 0.1 a.u. is generally sufficient for resolving spectra, but
at times propagation of the wavefunction requires finer steps).

1.5

Tools for analyzing HHG
To understand HHG, we must be able to analyze the effects of several different pathways.

Consideration of HHG spectra or time-domain pulses independently can be useful (e.g., to identify
or measure attosecond pulses), but insight into the underlying processes based on this information
is often limited. There are several possible tools to enhance our understanding of the mechanisms
and physics:(1) focus on a simpler approximate theory (e.g., strong-field approximation (SFA) [50–52],
classical trajectories), which can provide different access to what is going on than numericalsolutions of TDSE;(2) perform analysis in a mixture of the time and frequency domains; or,(3) widely scan laser and target parameters to seek trends, which can be especially useful given
how noisy HHG spectra often are.

20
When analyzing radiation emitted near the ionization threshold, many approximative methods fail due to the complex impact of excited states and atomic potential – which are, for example,
entirely ignored in the basic strong-field approximation [32, 50–52] and still poorly represented in its
extensions [53, 54]. Considering laser parameters that fall in the strong-field multiphoton (vs. thetunneling or perturbative multiphoton) regime of ionization and excitation similarly precludes theuse of most other approximate methods, or at least complicates them greatly, and therefore limits
intuitive insights based on the results (as in the case of generalized semi-classical trajectory methods, cf. Sec. 1.4.1). Therefore, we primarily utilize numeric solutions to the full time-dependent
Schrödinger equation (TDSE). The resultant spectra, however, often include many radiation peaks,
sidebands, and other features with sufficient magnitude to be of interest. On the other hand, thetemporal response during the pulse is consistently dominated by the low harmonics – i.e., one domain contains too much information and the other too little. Performing scans of laser (or target)
parameters will often at least clarify which radiation features are stable versus those existing only
in narrow slices of parameter space – and thus not likely to survive macroscopic propagation (seeCh. 4) or appear in experiment. This can illuminate the mechanisms behind a subset of featureswith simple dependence (e.g., FID signals which appear at the field-free |0i → |npi transitionenergies). Often, though, understanding mechanisms requires simultaneous parsing of time andfrequency information, and frequently in conjunction with parameter scans. In the following subsections, we will discuss the particulars of a few forms of time-frequency analysis and some remarks
on parameter scanning.

1.5.1

Temporal windowing

Perhaps the simplest way to gain insight into the interplay between time and frequency inradiation is to consider a windowed Fourier transform of the dipole acceleration,Z ∞P (ω) =−∞

2a(t)W (t)e−iωt dt ,(1.38)

21
with W (x) ∈ [0, 1] being some window function. We already introduced the commonly-used Blackman window WBlackman (x, T ) = 0.42 − 0.5 cos(2πx/T ) + 0.08 cos(4πx/T ) [55]. By convention, mostnumeric HHG spectra are reported using a Blackman window with width T equal to the laser pulse
duration. This eliminates radiation occurring after the end of the laser and generally smooths the
spectra.
However, windowing can also provide additional information, for example comparing thespectra windowed to the span of the pulse versus one including extended propagation past theend of the pulse. In a previous paper [3], we used this technique to demonstrate the impact of
populations of excited states with dipole-allowed transitions to the ground state on the spectra:
FID lines appearing as sharp peaks in the post-pulse radiation (see Fig. 1.7).
More generally, we may window to particular sections of the pulse (e.g., first half, peak ±2o.c., etc.) to identify the radiation emitted during different periods. The primary limitations of
windowing are loss of spectral precision – since fewer data points in the Fourier transform include
actual data – and identifying windows that provide useful information.

1.5.2

Spectral filtering and temporal profiles

In the other direction, we may consider the temporal profile of radiation for particular spectral
regions. For example, we could ask what pulse would be produced by filtering to the radiation near
the 9th harmonic. This is accomplished by Fourier transforming the dipole acceleration into thefrequency domain, applying a spectral filter, R, then inverse Fourier transforming back to the time
domain2
P (t; R) = F −1 R(ω 0 )F[a(t0 )] .

(1.39)

Unless otherwise noted, we use for such studies a filter with a two-harmonic flat top and223

intensity (arb. units)

1.0

(a)

0.8
0.6
0.4
0.2
0.0

7

intensity (arb. units)

89harmonic order (in 0.057 a.u.)
3

1.0

4 5 ...
(b)

0.8
0.6
0.4
0.2
0.0

7

8
9
harmonic order (in 0.057 a.u.)3

1.0
intensity (arb. units)

4 5 ...

4 5 ...
(c)

0.8
0.6
0.4
0.2
0.0

7

8
9
harmonic order (in 0.057 a.u.)

Figure 1.7: Radiation spectra generated during the pulse (red dashed lines) and those including line
emissions after the pulse (black solid lines) are shown for hydrogen driven by a 20 o.c. sine-squared,
800-nm pulse with peak intensities: (a) I0 = 3.4 × 1013 W/cm2 , (b) I0 = 6.0 × 1013 W/cm2 , and(c) I0 = 16.4 × 1013 W/cm2 . The vertical gray lines show the field free energy differences between
the np energy levels (up to 14p) and the 1s ground state. (Taken from [3].)23Gaussian falloff with σ = 0.1 harmonics, centered at harmonic N :

i
h


exp − 1 ( q−(N −1)) )2
q ≤N −1


2
0.1



RN (qω0 ) = 1
N −1≤q ≤N +1




h
i



exp − 1 ( q−(N +1)) )2
q ≥N +12(1.40)

0.1

Spectral filtering in this way can be quite powerful – especially when combined with parameterscans – but it is important to recognize artifacts of the filtering technique, including:• Fourier transforms and windows/filters involve a sometimes surprising number of factors
and phases which must be managed carefully to avoid unintentional shifts in time or magnitude scaling.
• If spectral filters are narrower than peaks or near-constant regions, it is possible to introduce
spurious frequencies, which can appear as temporal beating.
• Sidelobes resulting from too-sharp filter cut-offs, which may be appealing when attempting
to isolate narrow features.

1.5.3

Wavelet analysis

True time-frequency analysis provides information on the intensity of a signal in regions ofboth time and frequency, together. The simplest method is to scan the center of a temporal window
across the duration of data, then stack the windowed frequency response. This has a variety of
drawbacks including those already mentioned. Many are resolved using a more complex method,such as the continuous wavelet transform [56–58]. We specifically utilize a Morlet wavelet
max(t
Z E)

f (ω, t; tf ) ∝



2

ω
−iω(t0 −t)− 12 Ωω (t0 −t)
0ER (t0 )edt0Ω0(1.41)min(tE )where
Z ωmaxER (tE ) =0

ER (ω)eiωtE dω ,(1.42)

24
ω is the frequency of the wavelet, and t and Ω0 are the center and the variance of the Gaussian
window. The parameter Ω0 allows control of the time vs. frequency resolution of the wavelet. In
this thesis, we have chosen Ω0 = 6 or 54, which provide a good resolution in time or frequency,
respectively, for 800 nm lasers. We note that Eq. (1.41) provides the high harmonic spectrum forΩ0 → ∞:

Ztf
fHHG (ω; tf ) ∝0

a(t0 )e−iωt dt0 .

(1.43)−∞

While wavelet plots provide a great deal of information, there are drawbacks. Firstly, thedifferences in scale (e.g., between different harmonics) require careful selection of color axis limits
to avoid either hiding important features or enhancing noise to the point of obscuring details. Thisin turn leads to difficulties in comparing wavelet plots for different laser parameters. Even when
consistent and informative color scales are achieved, the sheer volume of information conveyed canmake comparisons visually complex.

1.5.4

Parameter scans

One of the strongest tools to understand any physical phenomena is to see what happens
under variation (“Does the heavier ball fall faster?”). In the context of theoretical physics, this isparticularly easy, as we may set any arbitrary laser or target parameters and analyze the simulation
results. In the simplest case of scalar observables, results may be plotted as a line against parameters(e.g., probability of ionization vs. driving intensity). For radiation, though, this is generally notas productive. For instance, plotting the intensity of a given harmonic obscures important features(red/blue-shift, sidebands, width) and becomes unstable to deviations from the expected behavior(e.g., when the “9th harmonic” is not centered at 9ω, or when nearby off-harmonic peaks havegreater intensity). For this reason, heatmaps are popular – allowing for plotting the full harmonic
spectrum as a function of some parameter. Similarly, we find great benefit in plotting scans of
the temporal profile of a given harmonic across, e.g., intensity. Finally, there are cases where theresults from a single simulation are already multidimensional, and scans can often not be visually25
processed without separate images (incl. slideshows and animated media). While important for
our research process, we avoid including such visualizations in this thesis.

1.6

Outline of this thesis
The rest of this thesis comprises three projects on high-order harmonic generation tied to-gether by an emphasis on investigation of near-threshold harmonics and their connection to atomic
characteristics (including the potential and excited states). Additionally, methods are presented
for efficient simulation and analysis applicable to broad parameter regimes.
In Chapter 2, we will discuss single-active-electron (SAE) potentials for atoms, ions, andmolecules. These potentials allow the use of existing single-electron methods (especially numeric
TDSE libraries) for simulation of more complex targets. We outline a method for generating SAE
potentials based on analytic forms fit to density functional theory potentials – without any parameters modified by hand. Using one such potential, novel agreement is shown with time-dependent
density functional theory calculations for simulation of a strong-field process, namely HHG. We introduce additional applications to HHG from multi-orbital atoms, and nonlinear polarization andionization of molecules.
In Chapter 3, we analyze the impact of field-shifted resonances on near-threshold HHG. We
focus on a particular set of spectro-temporal features associated with the 3p state of the hydrogen
atom shifting into 9-photon resonance, investigating the dependence on a variety of laser and target
parameters. We show that aspects of the behavior are quite general, and similar features appear
in a broad range of parameters, including more intense near-IR driving lasers, weak mid-IR pulses,
and helium driven by intense 400 nm pulses.
In Chapter 4, we introduce an efficient method of macroscopic summation built upon a manageable set of TDSE simulations. Specifically, the response of macroscopic gas jets are determinedthrough Monte Carlo summation of relatively few (105 − 107 ) point-like radiators. The acceleration
response from each radiator is interpolated from a small set (100-1000) of full TDSE calculations.
This method allows for accurate simulation of near-threshold HHG – for the thin medium or low

26
gas density regimes that are free from longitudinal phase-matching effects – which is not possible
with approximate methods such as the Strong-Field Approximation. The method is immediately
extensible to long wavelengths and to other atomic and molecular targets, and results are presented
for mid-IR pulses and helium driven by 400 nm pulses. We additionally demonstrate a distinct
interference pattern in the angular distribution of HHG for focal phase distributions associatedwith highly-focused pulses.
Finally, in Chapter 5, we conclude with a summary of the topics discussed within this thesis
and consider future directions of the field supported by this work.

Chapter 2
Single-Active-Electron Potentials1

A widely-used method to study nonlinear effects in larger laser-driven atoms involves numerical solutions of the time-dependent Schrödinger equation (TDSE) within the single-active-electron(SAE) approximation. In the SAE approximation, it is assumed that the laser field only interacts
with one – typically the most weakly-bound – electron in the atom or molecule, considering all otherelectrons as frozen spectators to the evolving physics. SAE-TDSE calculations have the advantage
of being a straightforward extension from common hydrogen TDSE codes. These calculations arecomputationally efficient and scale well to parameter regimes (e.g., long wavelength and/or highintensity) that are often infeasible to explore with multielectron codes. In cases where multielectron
calculations are possible, single-electron approximations still have value, since SAE-TDSE simulations can serve in initial exploratory parameter sweeps before running more expensive calculations
to identify multielectron effects.
Furthermore, ab-initio SAE calculations can provide important benchmark results for moresophisticated multielectron methods in parameter regimes where electron correlation effects are
small. In the perturbative interaction regime, such a comparison between results of time-dependent
density functional theory (TDDFT) simulations and ab-initio calculations using lowest-order perturbation theory has been done recently and a good agreement has been found [59]. This established
the applicability of TDDFT calculations for perturbative multiphoton processes. We are not awareof a similar study at higher intensities for highly nonlinear and nonperturbative processes. A few
1

Some results of this chapter are also presented in [1]28rare studies reporting a comparison between TDDFT and SAE-TDSE results for nonperturbativestrong-field processes typically attribute differences to the impact of multielectron effects [60, 61]
without establishing a benchmark test between SAE-TDSE and TDDFT calculations. The absenceof such benchmark tests partially contributes to the misconception in parts of the strong-field community that TDDFT studies cannot provide quantitative predictions for highly nonlinear processes
involving the continuum. In this Chapter, we show – by carefully constructing an SAE potential –that indeed agreement between SAE-TDSE benchmark results and TDDFT results can be found.
Hence, we establish a method to check the accuracy of TDDFT calculations in view of multielectron
effects and disprove the misconception about the quantitative predictive power of TDDFT calculations. We further extend the SAE framework to the generation of simple, but accurate, potentials
for molecular dimers. Throughout, we place an emphasis on technical obstacles we encounteredand advice on how others may similarly overcome them in future work.
In Sec. 2.1 we present the method to obtain analytical forms of SAE potentials based on DFT
calculations (which is essential for the benchmark tests). We establish the accuracy of the SAE
potentials by presenting predictions for energies of the ground state and singly excited states ofvarious atoms within the single-active-electron approximation (Sec. 2.1.4). In Sec. 2.2, we compare
results for high harmonic generation (HHG) from helium within the single-active-electron approximation with those from TDDFT calculations based on the same DFT approach for a few highly
nonlinear cases. Next, we show as another application of the SAE potentials results for high-orderharmonic generation from outer and inner valence shells (Sec. 2.2.3). Finally, we discuss SAE potentials for molecules and an application to nonlinear polarization and ionization in O2 (Sec. 2.3).
We summarize our results in Sec. 2.4 along with a brief description of published work by others
based on the SAE potentials.

2.1

Constructing SAE potentials
In this section, we discuss the construction of SAE potentials for atomic systems, based on

fitting to an effective Kohn-Sham potential. First, we discuss the concept in constructing the SAE

29potential for the benchmark tests before we briefly reiterate the density functional theory formalism
used to obtain the Kohn-Sham potential. Then we present the approach used to approximate theKohn-Sham potential via analytical functions and discuss strategies to determine the coefficientsin the model potential. Finally, we present the methods used to perform the ab-initio SAE-TDSE
and TDDFT calculations.
The construction of a potential to model single-active-electron behavior has been pursued
through several different strategies (e.g., [62–69]). SAE models are frequently developed using
analytic functions with free parameters that are fit to accurately simulate experimental behavior
in the scenario of interest (often by matching the ionization potential for the single active valence
electron to the experimental value, since many strong field processes – such as tunneling ionization– strongly depend on it). Besides the long-range behavior of the SAE potential – which is essential
for the determination of the ionization behavior – electron dynamics driven by a strong laser pulse
also probe the short-range properties of the potential (e.g., excited states). Constructing an SAE
potential that reproduces well ground, inner, and excited state energies, is therefore useful forstudying laser-driven processes from the valence as well as inner shells. Furthermore, it provides areliable ground for benchmarking TDDFT calculations against ab-inito SAE-TDSE results, as we
will demonstrate in Sec. 2.2.2. Conceptually, the SAE potential should therefore be based on thesame DFT approach used in the TDDFT calculations without any free parameter adjustment.
To this end, we have obtained analytical forms of single-active-electron potentials, using density functional theory (DFT) and physical reasoning regarding the electronic structure of an atomic
system. Fits to the effective Kohn-Sham potential are made without using ad-hoc free parameters.
We establish the accuracy of the potentials by comparing ionization and excitation energies. Basing
the approaches on the same functional we are set up to make an accurate quantitative comparison
between SAE-TDSE and TDDFT calculations for HHG from helium at intensities in the nonperturbative regime. We further note that by providing analytical forms of the SAE potentials, theSAE-TDSE vs. TDDFT comparisons can be reproduced or extended in the future.

30
2.1.1

Density Functional Theory

In the Kohn-Sham density functional theory (DFT) formulation of multielectron systems
KS (r) (i = 1, . . . N ) involved in the ground state of the spin-polarized[70, 71], each orbital ψi,σN -electron system is determined by solving a set of Schrödinger-like eigenvalue equations given by



1
KS
KS KS
− ∇2 + Vef f,σ (r) ψi,σ(r) = Ei,σψi,σ (r),2(2.1)

Vef f,σ (r) = Vext (r) + VH (r) + Vxc,σ (r)

(2.2)

where

is the effective Kohn-Sham potential and σ denotes the spin of the electron. The terms in Eq. (2.2)
R ρ(r0 ) 0
represent the external Coulomb potential Vext (r) = − Zr , the Hartree potential VH (r) = |r−r
0 | drPaccounting for the electron-electron repulsion with the total electron density ρ(r) = σ ρσ (r) =2
P PNσKS (r) , and the exchange correlation term V
ρ
(r)
with
ρ
(r)
=
ψ
i,σ
i,σ
xc,σ (r).
i,σσi=1
An exact expression for Vxc is not known and several popular approximations to this quantity
exist. We make use of the strategy proposed by Krieger, Li, and Iafrate [72] and extended in [73].
We choose a sum of local (spin-)density approximation (LDA) exchange and correlation functionals:
LDA
LDA
vxc,σ
(r) = µLDAx,σ (ρ↑ (r), ρ↓ (r)) + µc,σ (ρ↑ (r), ρ↓ (r))

where µLDAx,σ (r) = −16π ρσ (r)3(2.3)

and with µLDAc,σ , as defined in Appendix C of [74], and apply aself-interaction correction [74]:LDA−SIC
LDA
vxc,i,σ(r) = vxc,σ(r) − µLDAx,σ (ρi (r), 0)Zρi,σ (r0 ) 0LDA− µc,σ (ρi (r), 0) −dr .
|r − r0 |(2.4)

This results in an orbital-dependent potential, but the Optimized Effective Potential procedurefrom [72] resolves this via
KLI
Vxc,σ
(r) =

Nσ
X
ρi,σ (r)i=1

ρσ (r)(vxc,i,σ (r) + Xi,σ ) ,

(2.5)

31
where Xi,σ is the solution to the set of linear equations:S(δji − (Mσ )ji )Xi,σ = V̄xc,j,σ− v̄xc,j,σ ,

(2.6)

with
Z
(Mσ )ji =SV̄xc,j,σ=ρj,σ (r)Z
ρj,σ (r)ρi,σ (r)dr ,ρσ (r)NσX
ρi,σ (r)i=1

ρσ (r)(2.7)!
vxc,i,σ (r) dr ,

(2.8)

Z
v̄xc,i,σ =ρj,σ (r)vxc,j,σ (r)dr .

(2.9)

and note Xm,σ = 0 where m is the index of the highest occupied orbital (or indices for a degenKLI (r) has been found self-consistently byerate set), so those terms are ignored. The potential Vxc,σguessing a set of orbitals, i.e., ρi,σ , then using the orbitals to calculate an exchange correlation
potential (Eq. (2.5)), which is then used to update the orbital set as the lowest energy eigenstatesof the effective Kohn-Sham potential (Eq. (2.2)), repeating until convergence is reached. The DFT
calculations have been completed on a 1D spherical grid with spacing dr = 0.0025 a.u. and length
KS converges to a difference below 10−10
L = 80 a.u., requiring that the valence orbital energy Em,σa.u. between iterations. Note this DFT implementation was coded entirely independently of theOctopus [75–77] and Libxc [78] libraries used for our TDDFT comparison.

2.1.2

Notes on implementing DFTWhile we are by no means experts on density functional theory or its implementations, we
do have some observations which may be of use to other strong-field researchers venturing into
DFT. As with much of programming, the best solution is not to do it yourself. There are plenty ofpre-existing libraries (such as Octopus) that should suffice for most purposes. These are the resultof years of careful thought and optimization. Here, we chose to write a DFT code from scratch so
as to allow fully independent comparison with Octopus TDDFT. The following few tips may be
helpful when coding a new DFT implementation:

32• Many simple DFT algorithms function through iterative convergence. This necessitates
both that individual iterations have high accuracy, and that non-convergent behavior isavoided. An advantage, however, is that individual iterations can be manually adjusted,
as long as final convergence is achieved.
• As a simple example, any initial density distribution is valid. It may not converge, but it
will not converge to an incorrect answer.
• Especially for weakly-bound valence electrons (e.g., in lithium or negative ions), the density
tail of the DFT potential may extend further than expected. If the tail is not properly
constrained, small errors may compound, resulting in poor long-range behavior of the DFT
potential. If, for example, the density is too high far from the atom in one iteration, thepotential for the next iteration will be too shallow, resulting in even greater, incorrect
density. We managed this by using grids with a large extent and fixing the long range tailto an exponential decay for the initial and middle iterations. Additionally, the iterative
process was cancelled when the tail was found to be flat or growing rather than decaying.
• For our process, an exceedingly common result involved two- or three-iteration loops. Instead of converging toward a single energy, the energy would consistently switch between
two or three energies (and, presumably, corresponding states and potentials). Generally,none of these energies would be close to the correct one (estimated by other DFT results
or experiment, or compared after the fact). One method to avoid this involves mixing multiple previous iterations to set the initial conditions for the next one [79]. We did this at
the level of the density, mixing according to ρinput,i = (1 − α)ρoutput,i−1 + αρinput,i−1 . The
mixing parameter α limits the rate of change between iterations while also acting simply asa non-physical dial to seek convergent behavior – according to our experience, calculations
with α = 0.7 often converged and with decent speed, but it was frequently necessary toscan α to achieve convergence.

33• Computing a set of lowest eigenstates of some potential is a key component of the iteration.
The popular Imaginary Time Propagation (ITP) method [80] is ill-suited to this for anumber of reasons. Firstly, it is relatively slow, especially when computing eigenstates
above the ground state. Similarly, the accuracy falls off as higher eigenstates require manual
orthogonalization. Replacing ITP with the Arnoldi method [42] reduced the iteration timeby orders of magnitude, allowing for the generation of DFT potentials for atoms as large
as xenon. (Note: Arnoldi has similar benefits for generating initial TDSE states.)• Whether using a DFT potential directly or fitting to an analytic form, it is importantthat the short range behavior (generally a Coulomb singularity, with a constant offset as
seen in Figs. 2.1 and 2.2) is accurately captured. We made use of (constant) grids muchfiner than practical for time-dependent calculations and visually checked that Vxc,σ (r → 0)
approached a constant.

2.1.3

Analytic fitting

One can proceed by using the DFT potential numerically; it is however instructive and usefulto provide an analytical form of the SAE potential as:
VSAE (r) = Vlong + Vshort + Vshell ,(2.10)

where Vlong is the long-range Coulomb term:
Vlong (r) = −C0,r(2.11)and Vshort is a screened Coulomb (or “Yukawa”) potential
Vshort (r) = −Zc e−cr.
r(2.12)In test calculations we left C0 and Zc as free parameters in the fitting procedure. The results,
shown in Table 2.1, are close to the expected values of C0 = Z − (Ne − 1), where Z is the charge ofthe nucleus and Ne the number of electrons in the atom or ion, and Zc = Z − C0 accounts for theSpecies
He
Ne
Ar

C0
0.999441
1.066567
0.985250

Zc
1.000562
8.933208
17.014712

C0 + Z c
2.000002
9.99977717.999962

34

Table 2.1: Fitting results for charge parameters C0 and Zc in SAE potentials for rare gas atoms
helium, neon, and argon.

remainder of the charge. For this reason we did not vary these parameters in the remainder of ourstudies.
Vshell captures the effective SAE potential at intermediate distances with one exponential
term for each of the n orbital shells:
Vshell (r) = −nXai e−bi r .

(2.13)

i=1

This is motivated by the comparison of electron density and exchange correlation potential basedon DFT calculations (see Figs. 2.1 and 2.2) where we identify one smooth region per shell. We haveapplied the above potential form to both atoms and ions up to argon. For a given species the samepotential is found for both valence and inner orbitals. We note that for larger atoms, the present
approach based on an orbital shell picture may need to be modified, since in such atoms subshelldensities are more strongly overlapping (see Sec. 2.2.4).
The effective single-active-electron potential VSAE :n

VSAE (r) = −C0 Zc e−cr X −bi r
−
−
ai e
r
r

(2.14)

i=1

depends on n linear parameters (ai ) and n + 1 nonlinear parameters (bi , c). We fit VSAE (r) to
Vef f,↑ (r) using the Levenberg-Marquardt (LM) least-squares curve fitting algorithm [81, 82] implemented in the LMFIT library [83]. The LM algorithm involves methods to avoid local minima,
resulting in nonlocal fitting based on initial parameter guesses. As an extra check, the parameters
bi in the exponentials have been scanned across a range of initial guesses, selecting as final fit theone with lowest chi-squared value. In Table 2.2 we show the results for the parameters of our fits
for a number of atomic and ionic species. The relative error between the fit and the actual DFT

35

Figure 2.1: Results of DFT calculations for (a) electron density and (b) exchange correlation
potential of the helium atom (Taken from [1]).

Figure 2.2: Same as Fig. 2.1 but for the argon atom (Taken from [1]).

Zc
1
2
3
9
10
11
17
16
15
14
13
12
1
3
5
6
7
8
8
9
17
11
17
19
35
36
37
47
53

c
2.032915.9594
2.04810.88701.49271.42480.81030.86980.8792
0.94450.85290.89293.3926
1.98491.2052
1.37201.60001.78682.14161.07540.737516.00360.88268.92210.83561.08601.10390.72020.8608a1
0.39539.11240.5294-9.9286-11.3552-14.5892-15.9583-16.0391-16.4007-16.4800-17.4441-17.5407
0.76690.3089-4.4396-5.2033-5.4268-5.9589-6.1507-11.6217-14.1508
23.4319-17.8493
7.4822-30.9586-36.3604-39.7948
0.07760.10391.60281.81412.17682.17802.16812.15502.11352.08181.94573.73372.410439.6268-102.2050-122.5379-97.8013-61.0621-67.7122

1.84010.55160.69210.77370.90841.0293-5.8563-26.3671
98.1002-29.0486
56.3383-65.9241-58.1493-59.0982-32.4798-44.7989a3

19.31450.3219-5.9950-2.0302-1.9433-27.7467-26.9860-26.6805-25.8243-26.0893-25.4398

a2

2.038072.781394.571769.8723-95.5870-104.9326

a4

0*
0*
7.94449.8085a5

b1
6.18053.60400.84751.37462.55972.70011.23051.3146
1.34861.45211.41411.50248.39810.39102.38382.75813.1481
3.48883.83711.64531.13594.1702
1.33651.27981.20651.64391.73370*0*1.433548.6696
50.920460.7533
63.890769.9272
4.1576
4.0661
11.2679
4.7361
3.9896
3.3382
3.85304.03631.01781.214811.3082
37.5567
3.796310.1466
12.3150
4.39464.45144.46564.51714.46134.4823b2

96.3853
62.328177.55658.934718.6001
20.2399
20.4979
2.67483.404147.954251.7100
86.7179
88.331590.3068
94.5151101.5018
108.4695b3

144.0899
23.932224.7360
26.7700
10.7140
12.4774b4

0*0*148.2401
133.1390

b5

Table 2.2: Parameters for VSAE (Eq. (2.14)) obtained in analytical fits of single-active-electron potentials for various atoms and ions.
(The upper part of this Table also appears in [1].) (*) For Rb, Sr, Cd, and Xe, see Sec. 2.2.4.

He
Li
Be
Ne
Na
Mg
Ar
Ar+
Ar2+
Ar3+
Ar4+
Ar5+
Li+
Li−
C
N
O
F
Ne+
Na+
Cl−
Ar6+
K+
Ca
Kr
Rb
Sr
Cd
Xe

C0
1
1
1
1
1
1
1
2
3
4
5
6
2
0
1
1
1
1
2
2
0
7
2
1
1
1
1
1
1

36

37

Figure 2.3: Relative error between results of the analytical fit and DFT calculations for the
single-active-electron potentials for (a) helium and (b) argon (Taken from [1]).

Species
He
Li
Be
Ne
Na
Mg
Ar
Ar+
Ar2+
Ar3+
Ar4+
Ar5+

Experiment
0.9037
0.1981
0.34250.79270.18880.28110.5792
1.0154
1.4972
2.19802.75703.3446DFT
0.94530.19700.32620.83530.21510.29790.6003
1.07151.58502.15892.78353.4419

SAE
0.944
0.200
0.326
0.834
0.228
0.302
0.607
1.078
1.593
2.166
2.791
3.446

38

Table 2.3: Comparison of first ionization energies: Experimental data [2], DFT calculations, andanalytical potential fit (SAE) (Taken from [1]).

results for Vef f,↑ is presented in Fig. 2.3 for helium and argon.

2.1.4

Ionization and excitation energies

We have obtained energies for the ground state and (singly) excited states of the atoms andions for the analytical SAE potential fits using the Implicitly Restarted Arnoldi Method [43]. For
the analytical SAE potentials we have performed calculations on the same 1D spherical grid as forthe DFT calculations, with dr = 0.0025 a.u. and L = 80 a.u..
In Table 2.3 we compare the ground state energies for the atoms and ions studied, obtainedin the DFT calculations and with the analytical SAE potential fits, with the experimental datafrom the NIST Atomic Spectra Database [2]. The SAE potential results reproduce the DFT results
remarkably well, showing the success of the fitting procedure. The same conclusion holds for theionization energies of inner orbitals, as the example for neon atom in Table 2.4 shows. Further improvement concerning the agreement between the theoretical results and experimental data requires
more advanced DFT calculations (investigations suggest the next step up the functional ladder –the Generalized Gradient Approximation – is not sufficient to noticeably improve agreement). In
Tables 2.5 and 2.6 we present a comparison between experimental data [2] and the predictions

39

State
1s
2s

Experiment35.1590
1.8134DFT
30.8598
1.6723SAE
30.9051
1.6770

Table 2.4: Same as Table 2.3 but for ionization energies of inner orbitals in neon (Taken from [1]).

State
2s
2p
3s
3p
3d

Experiment-0.14595-0.12382-0.06126-0.05514-0.05561

SAE-0.15969-0.12847-0.064999-0.056679-0.055581

Table 2.5: Comparison of experimental data [2] and predictions using the analytical potential fit(SAE) for ionization energies of excited states in helium (Taken from [1]).

State
4s
4p
3d
5s

Experiment-0.14834-0.08937-0.05606-0.05581

SAE-0.17196-0.10843-0.07170-0.06463◦ )X argon excited states (Taken from [1]).
Table 2.6: Same as Table 2.5 but for 3p5 (2 P1/2

40
obtained using the SAE potentials for some singly excited state energies of the helium and argon
atoms, which shows a reasonable agreement.

2.2

HHG with SAE potentials

2.2.1

SAE-TDSE and TDDFT methods

To obtain HHG spectra in the SAE-TDSE calculations, we have propagated the wavefunction
using the static potential VSAE (r) and the time-dependent laser interaction potential described inlength gauge as:
VI (t) = r · E(t),(2.15)

where E is the electric field of the laser. In the present calculations, the initial state was set as
the ground state of the helium SAE potential. We have discretized the wavefunction and (radial)potential on a cylindrical 2D grid with dρ = dz = 0.03 a.u, and grid sizes up to ρmax = 50 a.u. andzmax = ±100 a.u.. We have applied the 2nd order Crank-Nicolson method and the split-operator
approximation to propagate the wavefunction starting from the initial state.
To study the accuracy of the SAE potentials, we compare results of SAE-TDSE calculations
with those of TDDFT calculations (performed by Tennesse Joyce). For the latter, the wavefunction
has been propagated using the Kohn-Sham orbitals obtained by solving the time-dependent KohnSham equations



1
∂ KS
KS
− ∇2 + Vef f,σ (r, t) ψi,σ(r, t) = i ψi,σ(r, t),
2∂t(2.16)

Vef f,σ (r, t) = Vext (r, t) + VH (r, t) + Vxc,σ (r, t) ,(2.17)

with

where
Vext (r, t) = −Z+ VI (t)r(2.18)and the Hartree and exchange correlation terms depend now on the time-dependent density ρ(r, t).
The TDDFT calculations have been performed using the Octopus code [75–77] on a Cartesian 3D

41
grid with spacing dx = 0.3 a.u. in a cylindrical box of length at least 70 a.u. in the z-direction, anddiameter at least 40 a.u. in the xy-plane. An additional 20 a.u. in every direction outside the box
was added for the complex absorbing potential. Note that use of 4th-order finite difference allows
for coarser dx here.
We have taken the laser to be linearly polarized in the z-direction, with a sin2 pulse envelope,variable central frequency ω0 , and full-width duration T = 8π/ω0 corresponding to four optical
cycles. That is,2

E(t) = E0 sin



πt
T


sin(ω0 t)ẑ

(2.19)

for 0 ≤ t ≤ T . We have evaluated the dipole acceleration a(t) using the Ehrenfest theorema(t) = h∇V i + E(t).

(2.20)

where V = VSAE (for SAE-TDSE) and V = −Z/r (for TDDFT). In the case of TDDFT, thedipole accelerations from all orbitals are added together, and other interaction terms cancel. The
harmonic spectrum P (ω) is then obtained by taking the windowed Fourier transform of the dipole
acceleration,
Z T
P (ω) =

2−iωta(t)W (t/T )edt ,(2.21)

0where W (x) = 0.42 − 0.5 cos(2πx) + 0.08 cos(4πx) is the Blackman window [55].

2.2.2

Comparison of SAE-TDSE and TDDFT results

The good agreement between analytical fit and DFT results provides the grounds to proceed with the quantitative comparison between TDDFT and SAE-TDSE calculations for HHG. In
Fig. 2.4, we compare results of SAE-TDSE and TDDFT calculations for the HHG spectrum driven
by a 4-cycle laser pulse with peak intensity of Ipeak = 1015 W/cm2 at central wavelengths of (a)λ = 267 nm (ω267 = 4.65 eV), (b) λ = 400 nm (ω400 = 3.10 eV), and (c) λ = 600 nm (ω600 = 2.07eV). For each comparison, we have matched the two results at one point, namely the energy of thedriving laser pulse. We used the helium atom as a target, this being the simplest atomic systemwith electron correlation effects. For laser wavelengths in the ultraviolet and optical regime, we42Figure 2.4: Comparison of TDDFT results (dashed red line) and SAE-TDSE (solid black line)
results for HHG spectra generated in the interaction of a helium atom with a 4-cycle laser pulse atpeak intensity of 1015 W/cm2 and central wavelengths of (a) 267 nm, (b) 400 nm, and (c) 600 nm.
Vertical dashed lines denote the ionization potential IP and expected HHG cut-off energy for eachcalculation. (Taken from [1].)43
expect that the correlation effects on the nonlinear and nonperturbative dynamics may be small.
Thus, the single-active-electron approximation should be applicable and, hence, constitute a test
bed for the comparison.
At the shortest wavelength (panel (a)), it requires six photons to exceed the field-free ionization potential Ip = 25.7 eV (cf. Table 2.3), while the energy of seven photons is needed to exceed
the field-shifted ionization threshold of Ip + Up = 32.3 eV where Up = 9.3 × 10−5 Ipeak [PW/cm2 ] λ2[nm] is the ponderomotive energy of a free electron in a laser field at the peak intensity. Furthermore, the HHG cut-off energy of 3.17Up + Ip = 46.7 eV ≈ 10 ω267 . The results of the TDDFT and
SAE-TDSE calculations are in excellent agreement over the whole spectrum, up to the expected
cut-off and even beyond. The position of the maxima and minima, the form of the odd high-orderharmonic lines, as well as the cut-off in the spectrum, are nearly identical in the results of both
independent calculations. We note small deviations between the results in the region of the 3rd and5th harmonics, which either may be related to differences in the excited state spectra for the two
calculations or may indicate minor electron correlation effects. However, overall the comparison
shows that electron correlation effects appear to be negligible for this nonlinear process in helium
at a UV driver laser wavelength.
Similarly good agreement is found at the longer central wavelengths of 400 nm (b) and 600
nm (c); the other laser parameters are held the same as for the calculations at 267 nm. We notethat while the 267 nm calculation was in the so-called multiphoton regime, the Keldysh parameterγ = 0.62 [50] indicates that the process at 600 nm occurs in the tunneling interaction regime. Dueto the longer wavelength, the ponderomotive energies are significantly larger. This results in larger
field-shifted ionization potentials and corresponding HHG cutoffs of 73.1 eV (b) and 132.2 eV (c),corresponding to the 24th and 64th harmonic of the driving frequencies, respectively. While thedifferences between the results of the two independent calculations are slighter larger than at theUV driver wavelength, they remain overall clearly less than a factor of two. The largest deviations
are again found in the energy region near the ionization threshold. We note that at higher energies,
even subtle features are in excellent agreement for the two calculations.

44
The agreement between the results demonstrates the accuracy and applicability of the singleatom potentials and shows the capability of TDDFT to provide reliable results for nonperturbativestrong-field processes in this parameter regime. The deviations in the ionization threshold region
may be attributed to a minor role of electron correlation and/or excited states during the process.
We note that studies reporting a comparison between TDDFT and SAE-TDSE results for nonperturbative strong-field processes are rare [60, 61] and we are not aware of another study in thestrong-field and tunneling regime reporting agreement in parameter regimes of small correlation
effects. Typically, differences between results have been interpreted as the effect of correlation
effects [60], even in parameter regimes in which we find excellent agreement between SAE-TDSE
and TDDFT results. This demonstrates the importance of the method in constructing SAE potentials used in the present study and the general benchmark tests performed. It also clarifies thesometimes expressed misconception that TDDFT calculation cannot provide accurate results in thehighly nonlinear multiphoton and tunneling regimes.

2.2.2.1

Subtleties of convergence

Informal discussions suggest that other groups have sought to compare TDSE and TDDFT
simulations of HHG and generally found poor agreement. We therefore think it is valuable toaddress the process behind our final results, including various obstacles.
First off, TDDFT calculations are generally more expensive than TDSE calculations. Obviously, the extra computations required to periodically update the DFT potential contribute; but
additionally, TDDFT codes usually have different goals than TDSE codes and are optimized correspondingly. In the case of linear polarization and atomic targets, most modern TDSE codes use
a spherical harmonic basis (as in our current TDSE code) or at least exploit azimuthal symmetry(as the code used for the comparisons in this Chapter) for speed. Since most targets in TDDFT
codes do not have spherical symmetry, the option is rarely included. Thus, even single-electron
calculations (e.g., hydrogen) will often be much more expensive when run with a TDDFT library,even though the actual physics is identical.

45
As a result of this expense – and the differing needs of common calculations – TDDFT simulations are often run with larger grid spacing and smaller box size than is common in TDSE
simulations. Consequently, HHG convergence tests may (and, in our case, did) begin with the
TDDFT parameters far from being converged. This led to more complicated interactions between
simulation parameters than is typical for almost-converged TDSE cases. As one example, increasing the box size along the laser polarization led to worse agreement! After further investigation,
we came to the conclusion that the issue was related to the relatively narrow box perpendicular topolarization – electron wavepackets were artificially reflecting off the boundaries, leading to unexpected interference behavior as a function of box length. After increasing the box width sufficiently,
changes in the box length displayed a more normal convergence behavior.
On a different note, HHG depends significantly on the ionization potential. When this does
not agree between SAE-TDSE and TDDFT simulations, the spectra will differ. This suggests first
the utility of basing the two simulations on a shared potential formulation (e.g., a particular DFT
functional, in our case with one being static and fit to an analytic form). As above, it is also
important that both codes be converged for the ground state energy (here, grid spacing is mostimportant). Some of these notes may seem obvious, but when the result is a month-long TDDFT
calculation, there is the incentive to cut corners and assume that the simulation is sufficiently
converged when it may not be.
Finally, we provide a few notes relevant to any comparison of HHG calculations, including
those using different TDSE implementations. In principal, HHG can be computed equivalently via
the dipole moment, the dipole velocity, or the dipole acceleration. Since the dipole moment is thesimplest to compute, this approach may look most attractive. However, the equal treatment ofelectron density anywhere on the grid leads to errors when ionized wavepackets must be absorbed.
Using the dipole acceleration (specifically via the Ehrenfest theorem, as shown above) minimizes
these issues and we have generally found that method to provide better convergence behavior withrespect to various numeric parameters, especially box size. The quality of the HHG spectrum using
the dipole velocity falls between the other two, but we have generally found its accuracy to be46
closer to calculations based on the dipole moment. Thus, we highly recommend the use of dipole
acceleration, even if adding it to a library appears complicated. For a thorough investigation anddiscussion of more theoretical background, see [84].
Another potential issue can arise when a team is using multiple libraries since it is easy for
different definitions of parameters to come into conflict. Often, the results will only differ slightly,
making it feasible that the difference is intrinsic to methods rather than a result of miscommunication. A few examples:• Is the “length of the box” defined by diameter or radius?
• How, exactly are the laser wavelength and frequency defined? For example, is the wavelength exactly 800.0 nm or is the frequency set to 0.0570 a.u.?
• Is the duration of the pulse measured in FWHM (Gaussian) or full width (sine-squared)?
Are reported durations in fs rounded or does the pulse have a non-integer number of cycles?
• Is the laser field defined by an expression for electric field or for vector potential? The
difference introduces a CEP offset and frequency shift [85], both especially relevant for
short pulses. This is independent of gauge as the other field may be derived throughnumeric or analytic differentiation/integration.
Lastly, one needs to be careful to ensure consistency of HHG calculations using the dipole
acceleration. Perhaps the simplest method is to pass the raw acceleration data through the sameplotting script. Otherwise, one needs to be aware of the following considerations:• Is the acceleration windowed before applying the Fourier transform? If so, what window is
used?
• Are the Fourier transform conventions consistent? If not, an overall factor may appear.
• Do signals have consistent time sampling and padding? This is mostly a technical issue forplotting, but important to consider when things go wrong.

47
2.2.3

Multiple SAE-TDSE for larger atoms

Beyond the application for benchmark-testing TDDFT calculations, the present SAE potential can of course be applied for other studies as well. Recently, the influence of inner valenceshells on high harmonic generation has been studied [86–89]. The effects in the spectra are often
interpreted by multielectron or intershell correlation, so accurate single-active-electron calculations
can serve as a useful tool for comparison and, hence, identification of the correlation effects.
Using the potentials for neon, we have obtained high harmonic spectra from the outermost
2p shell (m = 0) and the inner 2s valence shell. In the calculations, we have applied a laser pulseof 4 optical cycles at 800 nm and a peak intensity of 1014 W/cm2 . Both spectra (cf. Fig. 2.5(a)and (b)) show the characteristic features of an HHG spectrum with plateau and cut-off. However,in both spectra there occurs a prominent peak between the 15th and 17th harmonic, at the energy
corresponding to the energy difference between the 2p and 2s states in our SAE potentials for neon.
The peak is an artifact of the SAE-TDSE calculations which allow for transitions between the subshells, although they are filled in the real multielectron atom. However, the unphysical feature iscompletely removed once the amplitudes for the processes driven from the two sub-shells are added
and the full HHG spectrum (Fig. 2.5(c)) is obtained. This is another indication of the accuracyof the single-active-electron potentials obtained with the present method, which therefore can beapplied for single-active-electron reference calculations, e.g., for the identification of multielectron
correlation effects.

2.2.4

Notes on SAE for large atoms

For larger atoms, the above analytic form, Eq. (2.14), does not always capture the shape of
the DFT potential well. The resulting poor quality fits are generally easy to spot from the outputof the fitting function. A few red flags are:
• The error is much higher than for similar atoms.
• Parameter values fall very close to the limits of the parameter range.

48

Figure 2.5: Results of SAE calculations for high-order harmonic spectra generated from (a) 2p-shelland (b) 2s-shell in neon by a 4-cycle laser pulse at peak intensity of 1014 W/cm2 and wavelengthof 800 nm. The vertical dashed line indicates the 2s − 2p transition energy at which an artificial
peak occurs in the spectra from the two sub-shells, which is removed once the the amplitudes forboth processes are added (Panel (c)). (Taken from [1]).

49• Parameter values for different terms are close together, or nearly opposite. For example,
we have encountered cases where b1 ≈ b2 and a1 ≈ −a2 .
• Correlations between terms are reported as 1.0.
• The uncertainty in fitted parameter values is high (esp. when > 100%).
• Other warnings from the program, such as “Reached maximum iterations without achieving
requested tolerance.”
Often, these issues can be resolved by modifying the fitting options, such as increasing theallowed parameter ranges, the number of iterations, or the initial parameter guesses. But, whenproblems persist, it is useful to reconsider the analytic form.
For our SAE potentials, we have identified a few species where the predicted number of
exponential terms overfits the data. Specifically, neutral Rubidium and Strontium (with 5s valenceshell) produce poor fits with five exponential terms. Aside from a high fit error, the corresponding
fits display nearly every red flag mentioned above. However, when fit with four exponential terms,
the issues disappear while the fit error remains reasonable. These are the fits reported in Table 2.2.
For the two largest atoms considered, cadmium and xenon, we encountered different behavior.
Namely, one of the “exponential” terms in fact fit to a constant (i.e., b1 = 0.00 with nonzero a1 ).
Interestingly, no other red flags appeared. Reducing the number of terms to four or three decreased
the fit accuracy but did not eliminate the constant terms. Physically, constant potential terms
are not reasonable but it is unclear whether the issue is related to the DFT potential or that ouranalytic form is missing some characteristic of these large atoms.

2.2.5

Effects of manual SAE parameter modification

Since many previous SAE potentials have included manual modification of parameters tomatch experimental ionization potentials, we briefly explore the effects of various modifications of
VSAE on a helium HHG spectrum. For this particular investigation, we included only the exchange

50
terms in Eqs. (2.3-2.4):
LDA−SIC
LDA
vx,i,σ(r) = µLDAx,σ (ρ↑ (r), ρ↓ (r)) − µx,σ (ρi (r), 0) −Zρi,σ (r0 ) 0dr .
|r − r0 |(2.22)The resulting IP values are lower than those for the vxc functional by up to about 5%, and fitting
parameters also vary. However, we do not expect this to affect the manner in which HHG changes
as a function of parameter values.
To begin, we have investigated the effects of removing individual terms from VSAE andmodifying other parameters to keep IP constant. We have performed the matching on the actual
2D cylindrical grid used for HHG calculations, with radius R = 40 a.u., length L = 120 a.u., andgrid spacing dx = 0.1 a.u.. Due to limitations of convergence in this older code, IP = 0.9477 differsfrom the converged value for the eigenstate of the potential, −0.918 (consistent for numeric Vef f
and fitted VSAE , on the finer grid described in Sec. 2.1.4). We have determined the parametervalues necessary to match IP – to at least three decimal places – using a simple secant search.
In Table 2.7, we show six such parameter sets with parameters and eigenstate energies on thecylindrical grid. The first three potentials share nearly identical energies for all states, while theremaining three differ greatly (up to 140%), or lack bound excited states entirely.
In Fig. 2.6, we show HHG with a driving laser of 1015 W/cm2 , 4 o.c., and 800 nm for each ofthese potentials. The results are indistinguishable for the first three, so we merge the legend entries.
After the 80th harmonic, the spectra show decent agreement – despite significant alterations to:the short- and long-range behavior of the potentials, and the excited state spectra. In this highly
non-perturbative ( UωP = 38) and tunneling-dominated (γ = 0.47) regime, the value of the ionization
potential appears to be the sole target feature influencing the shape of higher harmonics. However,radiation near and below the ionization threshold is strongly impacted by the atomic potential(and, as discussed in Ch. 3, excited state spectrum). This is especially clear for the short-range
potential, YE-a, which lacks excited states. In this case below-threshold radiation is strongly
suppressed relative to the plateau harmonics, as well as compared to below-threshold radiation ofother potentials.

1
− −
r

1
− − a0 e−brr
102.3526
7.09104.0762-0.948-0.159-0.128-0.065-0.056-0.056-0.035-0.032-0.031-0.031r
1
12.3526
0.66284.0762-0.948-0.160-0.128-0.065-0.056-0.056-0.035-0.032-0.031-0.031− ae−br

e−cr
1
− −
r
r
1
1
2.064004.0762-0.948-0.160-0.128-0.065-0.057-0.056-0.035-0.032-0.031-0.031

0
e−c r

CY-c
1
0
− − ae−b rr
102.3526
0.66280.3130-0.948-0.354-0.306-0.153-0.131-0.113-0.057-0.057-0.053-0.049CE-b−0
12.352616.8190
4.0762-0.948
n/a
n/a
n/a
n/a
n/a
n/a
n/a
n/a
n/a

r

− a0 e−br

YE-a
e−cr

C00
r
1.368202.064004.0762-0.948-0.236-0.234-0.105-0.104-0.104-0.059-0.059-0.059-0.059−

C-C0

Table 2.7: Modifications of helium SAE potential keeping IP = 0.948. In the first row, potentials are labeled by which terms are included(Coulomb, Yukawa, Exponential) and which term is modified from the base parameters. The second row defines the form of the modified
potential, with specific parameter values given in the following rows. Finally, the eigenstate energies are given through E4f . The potential
without a Coulomb term (YE-a) does not support bound states above |1si. We note the consistency of state energies (and, as will beseen, HHG) for the first three potentials.

C0
Zc
c
a
b
1s
2s
2p
3s
3p
3d
4s
4p
4d
4f

VSAE

CE-a

CYE

51

52

Figure 2.6: Results of SAE calculations for high-order harmonic spectra generated in the interaction of a helium atom with a 4-cycle laser pulse at peak intensity of 1015 W/cm2 and central
wavelength of 800 nm. Each line represents results for an SAE potential defined in Table 2.7. The
black line corresponds to the original potential as well as two other potentials with indistinguishable
spectra. Vertical dashed lines denote the shared ionization potential (IP = 0.948 a.u.) and HHG
cut-off (8.0 a.u.). Results have been normalized to the intensity at the cut-off.

53
Additionally, we have performed calculations in which we have varied the SAE potential
to match the experimental ionization potential. These HHG spectra (not shown) differed fromthose in Fig. 2.6, but the areas of agreement between different potentials and qualitative differences
remained consistent. We note that the modified full potential produced indistinguishable stateenergies and HHG regardless of whether the a, b, or c parameter was chosen for modification.

2.3

Molecular SAE2
Here, we extend the SAE formulation to the oxygen molecule for investigation of nonlinear

polarization and ionization. In particular, we collaborated with the group of M. Kolesik (Universityof Arizona) in extending their previous applications of the Metastable Electronic State Approach(MESA) [90] to O2 . In a previous paper, their group investigated molecular nitrogen using multielectron calculations [91], but attempts to use a simplistic model potential for O2 did not achieve
the expected level of agreement with experimental data. Our construction of the molecular SAE
potential follows primarily the same procedure as that for atomic SAE potentials above.
The effective Kohn-Sham potential Vef f,σ (r) = Vext (r) + VH (r) + Vxc,σ (r) is here calculatedusing the Octopus code (our separate DFT implementation is limited to species with radial symmetry) on a cylindrical 2D (ρ, z) grid with spacing 0.05 a.u. and length of 10 a.u. in both dimensions.
The VH + Vxc,σ terms have the same form as in Sec. 2.1.1 above, now with two Coulomb terms forthe interaction with the nuclei, so
Vext (r) = −88−.
|r − r1 | |r − r2 |(2.23)As above, all core electrons are included in the calculations (i.e., no pseudopotential is used). Here,the valence electrons are spin-polarized.
For the analytic form, we use a pair of radially-symmetric functions centered at each nucleus
plus a cylindrically-symmetric term along the molecular axis. This approximately captures thedetails of the potential from the core and valence electrons while retaining a simple form. Thus,
2

The results of this Section are also presented in [5]54
for a molecule aligned with the o-axis, we fit an analytic potential of the form
V (o, p, q) = f

p

p

 p

(o − L/2)2 + p2 + q 2 + f(o + L/2)2 + p2 + q 2 + g o, p2 + q 2 , (2.24)

where L = 2.3 a.u. is the bond length,
f (r) = −C0 Zc e−r/c−+ a1 e−r/b12r2r(2.25)

with the net charge seen by the active electron, C0 = Z − (Ne − 1) = 1 split evenly across thetwo nuclei, as well as the shielded charge, Zc = Z − C0 = 15, where Z (Ne ) is the total number ofprotons (electrons) in the molecule, and22g(o, ρ) = de−(o/co ) −(ρ/cρ )(2.26)

with the parameters obtained from a least-squares fit, c = 1.326, a1 = 10.804, b1 = 0.466, co =3.248, cρ = 2.778, and d = 1.643 (all in Hartree atomic units).
The radial functions f (r) are a truncated form of the above atomic SAE potentials, dropping the second (n = 2) exponential term for efficient fitting. Note that the parameter definitions
are different from above (for one, nonlinear parameters are here given as distances). The additional cylindrical Gaussian term g(x, ρ) with three fitting parameters coarsely captures the effectof molecular bonding orbitals.
The numerical effective Kohn-Sham potential from the Octopus code is fitted to the analytic
form (2.24) using a similar least-squares fitting procedure to that above. Fitting of the radial
term is facilitated by (Fourier) transforming the 2D potential to a radial form. This improves
the speed of the fit, but is not strictly necessary. The accuracy of the fit can be seen from thecomparison of the analytical SAE fit with the numeric DFT potential in Fig. 2.7. Shown are cuts
along (a) and perpendicular to (b) the internuclear axis, as well as the relative error as a functionof the coordinates (c). Since the electron-nucleus potential Vext is the same in the two forms of
the potential, we have excluded it from the visualization to enhance the comparison. We further
note that the ionization potential of the analytic fit (Ip = 0.43 a.u.) agrees well with that of theOctopus calculations (IP = 0.40 a.u.).

55

Veff, σ − Vext
V − Vext

15
10

(a)
−10

−5

0

o (a.u.)

5

10

8
6
4
2
0

(b)
2

4

0.15

4

0.10

2

0.05

0

0.00

−2

−0.05

ρ (a.u.)6

ρ (a.u.)8

10

Relative Error

5

Veff, σ − Vext
V − Vext

10

potential (a.u.)potential (a.u.)

20

−0.10

−4 (c)−5.0 −2.5 0.0

2.5

o (a.u.)

5.0

−0.15

Figure 2.7: Top row: Comparison of the analytical SAE potential fit (dashed red line) and thenumeric DFT potential (solid black line) (a) along the molecular axis (ρ = 0) and (b) for a line
bisecting the molecular axis at the origin (o = 0). Bottom row: (c) Relative error between numeric
DFT potential and analytic SAE fit (Eq. (2.24)). Note that the error, within a range of ≈ 1.5 a.u.
from the nuclei, is ≤ 5% (Taken from [5]).

H2
O2

L
1.43142.3000

C0
1
1

Zc
1
15

c
0.316
1.326

co
13.8413.248

cρ
18.6372.778

d-0.162
1.643

a1
1.205
10.804

b1
0.497
0.466

56

Table 2.8: Parameters for molecular VSAE (Eq. (2.24-2.26)) obtained in analytical fits of singleactive-electron potentials for various molecules.

2.3.1

Potentials for other molecules

We have also begun to construct analytic fits to other molecules – thus far, H2 based on
Octopus DFT calculations completed by Lauren Bauerle. For this and later symmetric dimer
molecules, we use the same analytic form, Eq. (2.24-2.26). The fitting process is unchanged andwe achieve similar results. For future work on asymmetric dimers such as CO, we suggest several
modifications:(1) Consider additional or different bridge terms gj .
(2) Utilize separate parameters (excl. bond parameters L, d, co , cρ ) for each atom.
(3) Allow charge coefficients C0,i and Zc,i to vary in the fit, while retaining a fixed sum.

2.4

Summary
To summarize, we have presented analytical forms for atomic and molecular SAE potentials,

which are based on results of DFT calculations. The potential fits are shown to provide a good
agreement with the DFT results for the energy of the ground state as well as with experimental data
for the energies of ground and excited states, and inner shells, of different atoms and ions, without
any ad hoc adjustment of parameters in the analytical forms. Using the potentials, a benchmark
test between SAE-TDSE and TDDFT calculations has been performed and a remarkable agreement
throughout the whole HHG spectra at ultraviolet and visible wavelengths has been achieved. Since
the calculations are based on two different approaches and different numerical codes, the agreement
indicates the accuracy of the SAE potentials as well as the capability of TDDFT to provide reliable57
results for a nonperturbative, highly nonlinear process such as high harmonic generation in thisparameter regime. Additionally, results of SAE calculations for HHG spectra from different shells
in neon reveal an artificial large peak at the transition energy between the 2s and 2p-states, whichis removed once the amplitudes of the processes from the two sub-shells are added up. Finally, we
discuss construction of SAE potentials for a usage beyond HHG, in the nonlinear polarization andionization of O2 .

2.4.1

Usage in other work

In addition to projects with my direct involvement (above), other members of our groupsuccessfully use the generated potentials for a variety of investigations, including on EUV ionizationof helium [92] and IR ionization of alkali and alkaline earth metals [93]. With the recent publicationof parameters and demonstrated agreement in helium HHG, we contribute to an ongoing discourse
on HHG from more complex atoms and molecules, including the first external citation [94].

Chapter 3
Resonance-Enhanced High-order Harmonic Generation

In the near-threshold region of the spectrum, off-harmonic spectral features, red-shifted withrespect to the harmonics, have recently received new attention in experiment and theory [39, 95–101], since they cannot be interpreted using the standard three-step picture of high-order harmonic
generation. Instead, different interpretations based on resonant excitation followed by free induction
decay [39, 96, 99], resonant enhanced harmonic generation [100], as well as population of excitedstates through frustrated tunneling ionization followed by free induction decay [101], have been put
forward.
In this Chapter, we investigate this spectral region with a particular focus on a set of radiation
features we associate with enhancement by the field-shifted 9ω → 3p resonance in hydrogen atoms.
In Sec. 3.1, we introduce this feature through spectral results for an illustrative case, and discussionof variations within the parameter space. Sec. 3.2 summarizes the prior work and semi-classicaltrajectory hypothesis of another group. We then demonstrate a method for identifying and verifying
radiation associated with enhancement from a particular resonance and its associated delay(s) inSec. 3.3 before applying these insights to several distinct features, which appear across a wider
range of driving intensities (Sec. 3.4). Next, in Sec 3.5 we return to the trajectory hypothesis anddemonstrate that it explains some but not all observations. In Sec. 3.6, we investigate deviations
from our model by keeping our core feature (resonance time) constant through simultaneous scans
across multiple parameters. In the remaining Sections, we provide a brief discussion of results inthe distinct parameter regimes of: short and very long pulses (Sec. 3.7), the helium atom driven

59
by a 400 nm laser (Sec. 3.8.1), and hydrogen atoms driven by a mid-infrared pulse (Sec. 3.8.2).
Finally, we provide a summary and discussion of future prospects in Sec. 3.9.

3.1

Spectral features
In Fig. 3.1, we show an example of a spectrum with harmonic peaks at ω, 3ω, 5ω, 7ω, and11ω, in which the radiation near the 9th harmonic is strongly suppressed relative to both nearby
frequencies and neighboring harmonics. Additionally, there is an extra peak below the 7th harmonic(near the 1s → 2p transition energy of 6.58ω) with similar intensity. These radiation features occur
during the pulse, and do not have the narrow width characteristic of field-free decay from excitedstates (cf. Fig. 1.7).
In order to gain further insights we have performed a series of calculations, for different peakintensities, in which the wavefunction has been propagated during the pulse and over a period offour times the pulse duration after the pulse. The results are shown in Fig. 3.2 as a function ofpeak intensity. It is seen that the spectral structure below the 9th harmonic, observed in the results
shown above, is present over a certain range of intensities at about 3.8 × 1013 W/cm2 . However, thestructure can be clearly distinguished from the free induction decay lines arising due to populationof excited states after the pulse ends – in Fig. 3.2(b) denoted and marked by red dashed lines.
In Fig. 3.3, results from similar sets of calculations varying driving laser wavelength (toppanel) and pulse duration (bottom panel) reinforce this conclusion, while highlighting trends in theenergy of the off-harmonic peaks. As either parameter is decreased, the side peaks clearly separate
further from 9ω. For shorter wavelengths (than about 790 nm), the 9ω radiation is dominant, but
the sidebands are still present. For the shortest pulses (less than about 7 o.c.), other structuresappear more intensely than both on-harmonic radiation and the feature present at longer durations.
Comparing spectral intensity and phase as a function of intensity in Fig. 3.4, the off-harmonic
feature appears to have less phase variation than the weaker lower-energy features, and to be related
to the on-harmonic emission at higher intensities. Combined with its significant strength within amoderate volume of parameter space, one would expect the off-harmonic radiation to be visible in60

Figure 3.1: High harmonic spectrum generated in atomic hydrogen, exhibiting a rich structure in
the near-threshold region of the spectrum between harmonics 5 and 12. The hydrogen ionizationthreshold of 0.5 a.u. = 8.77ω, where ω = 0.057 a.u. corresponds to the central driver wavelengthof 800 nm. Other laser parameters: sin2 pulse of 20 optical cycles with peak intensity of 3.8 ×
1013 W/cm2 . Note that the intensity of the spectrum is scaled by the signal at the fundamental
frequency.

61

Figure 3.2: HHG spectra as a function of peak intensity for (top) radiation during the pulse (asin Fig. 3.1) and (bottom) for propagation of the wavefunction over a period of four times thepulse duration after the end of the pulse. Results for each intensity have been normalized to the
maximum signal in the region between 8.75ω and 9.5ω. Red dashed lines denote field-free excitedstate transition energies.

62

Figure 3.3: HHG spectra as (top) a function of the driving wavelength at Nτ = 20 o.c., and (bottom)
a function of pulse duration at λ = 800 nm. For both plots, driving intensity I0 = 3.8×1013 W/cm2 .
Results include spectral contributions (free induction decay) for propagation of the wavefunction
over a period of four times the pulse duration after the end of the pulse. Results for each wavelength(top) or each duration (bottom) have been normalized to the maximum signal in the region between
8.75ω and 9.5ω. Red dashed lines denote field-free excited state transition energies.

63
experiment (this is further supported by macroscopic calculations in Sec. 4.4.2).

3.2

Resonantly-initiated quantum trajectories
A similar structure in HHG from argon was previously explained by Camp et al. in Ref. [100]via trajectories initiated by a resonantly-enhanced response through Stark-shifted excited states.
In pulses of moderate duration (e.g., 20 o.c.), the envelope varies slowly enough to consider
adiabatically-shifting states, but rapidly enough that the envelope cannot be considered constantover the scale of the travel time of an HHG trajectory. Assuming that excited states shift with~2

UP (t) = |E(t)|as described in Sec. 1.3, the resonance condition for sin2 pulses is given by:
4ω 2


I
Tr π
4
cos ±
= Nγ − (En − E0 ) ,
4ω 2
Nτ

(3.1)

where Tr is the resonance time in o.c. relative to the envelope peak, Nγ is the number of photons
to be absorbed for the transition, En is the field-free energy of the resonant state, I is the peaklaser intensity, ω is the central frequency, and Nτ is the pulse duration (full width) in o.c.. Intheir analysis, the radiation occurred at a specific delay after the resonances, td , measured in
optical cycles, dependent on the particular resonance but not on laser parameters. This delay was
attributed at least partially to the time associated with a semi-classical trajectory.
Camp et al. then claim the slope of the laser envelope is imprinted through the phase αUP /ω,which depends on time of recombination – thus resulting in a red-shift (blue-shift) when recombination occurs on the falling (rising) side of the pulse:∆ω = αdUP /ω
dt
t=±Tr +td

(3.2)

The scale of the phase is left as a free parameter, α, used to fit the predicted shift to theobserved radiation. We note that two other parameters – the degree of state shift as a fraction of UP ,and delay between resonance and emission – are modified to achieve agreement between temporalprofiles of radiation and the predicted emission time. The values of these free parameters are
justified as falling roughly within the range of previous studies (esp. on semi-classical trajectories).

64

Figure 3.4: HHG spectral (top) intensity and (bottom) phase as a function of peak intensity forradiation during the pulse (other parameters as in Fig. 3.1).

65
Our own investigations, below, confirm the association of temporal radiation with delays
after field-induced resonances (without by-hand modification – setting ∆E = UP (t)). However, inSec. 3.5, we find features not explained by the ∆ω prediction and present results suggesting thatsemi-classical trajectories do not apply at the parameters of our primary attention (specifically,
intensities roughly at and below 3.8×1013 W/cm2 ). Interestingly, we do find evidence of trajectorylike radiation at higher intensities.

3.3

Identifying resonances with temporal profiles
By windowing the spectrum and inverse-Fourier transforming, we may analyze the temporalprofile of particular harmonics. Here, the dipole acceleration is first Blackman-windowed and thena flat-top spectral window of width 2ω with Gaussian fall-off of width 0.2ω is applied. We identify
resonance times as a function of laser parameters according to Eq. (3.1). In the Figures below, wedraw these as lines connecting the times for discrete values of a scan. This generally leads to gaps
near the point where rising- and falling-side resonances meet (e.g., when the exactly-resonant peakintensity condition is not sampled), but allows for marking resonances on arbitrary scans of one ormore parameters.
A central benefit of these visualizations is the ability to match radiation features to particularresonance features. This is done by comparing the timing of radiation to that of various resonances.
While several resonances may occur during a given pulse, we will generally be able to identify asingle specific resonance as only this one will consistently follow the timing of the radiation acrossall parameter scans (and usually just 1-2 scans are sufficient). We then denote these resonances as
Nγ ω → |ni – for example, 9ω → 3p. Additionally, we often include the delay between resonanceand radiation, td , in optical cycles and which side of the envelope peak the resonance occurs, inphrases such as: rising 9ω → 3p +1.8 o.c. or f alling 9ω → 3p +0.5 o.c.. For the radiation featuresin Figs. 3.1-3.3, we will now outline the steps of identifying the associated resonance (9ω → 3p) anddelays (0.5, 1.8, 3.0 o.c.). In later sections, we primarily present the results of this process, where
associated resonances and corresponding delays have already been identified.

66
To begin, it is useful to plot a temporal profile scan as a function of driving pulse durationNτ , as shown in Fig. 3.5. In these scans, most dimensionless parameters (UP /ω, IP /ω, HHG cut-off,etc.) remain constant, as does the set of resonances which occur during the pulse. The timing ofresonances Tr does change, as measured in optical cycles or absolute time (e.g., fs), though not in
fractions of the pulse duration. In Fig. 3.5, we identify four temporal features of the 9th harmonic.
This radiation is primarily generated through a feature on the trailing end of the pulse. A somewhat
weaker contribution occurs on the rising edge, which seems to merge with the strong feature in the
trailing part for short pulses. Additionally, for longer pulses (more than about 8 and 15 cycles,
respectively) there are two additional radiation features close to the pulse peak, varying in relative
intensity as a function of pulse duration (strongest near Nτ = 15 and 25 o.c., respectively). To theeye, the slopes of all these features in the Figure, except for the strong feature on the trailing edge,
appear to be similar.
Next, in Fig. 3.6(a), we indicate the resonance times Tr for different transitions via dashedlines in the same plot as shown in Fig. 3.5. According to Eq. (3.1), Tr is proportional to the cycle
duration Nτ . Therefore, resonances appear as straight lines with slopes depending on the resonancechosen (and the other, fixed, laser parameters). As such, it is fairly simple to check which resonancelines match the slope of radiation features and identify the delay td . In this case, the central six
lines – (9ω → 4p), (9ω → 3p), (8ω → 2s), rising and falling – seem promising as matches to theslopes of the various radiation features. Upon checking, it becomes obvious that 9ω → 3p agrees
quite well with all features, while the other resonances do not. At the same time, it is possible toestimate the delays between this resonance and the radiation, finding td ≈ 0.5, 1.8, 3.0 o.c..
In Fig. 3.6(b), we draw only the 9ω → 3p resonances, with delays td . Contrary to our
expectations, it appears the rising 9ω → 3p resonance produces radiation at three distinct delays.
While the difference between the consecutive delays (∆td = 1.3 and 1.2 o.c.) is approximately equal,
our investigations have uncovered no further examples to mark this as more than a coincidence.
(For example, we have not found radiation matching 9ω → 3p + 4.25 o.c., nor other resonances
associated with radiation at more than a single delay.) As a side note, the unexpected result of these

67
multiple delays led us through a detour investigating the effects of pulse envelope and radiation inlonger pulses, which will be discussed in Secs. 3.5.1 and 3.7, respectively.
After identifying the resonance and delays associated with the radiation features in the present
results, we verify our conclusions by scanning other parameters. Firstly in Fig. 3.7, we present results obtained by varying the driving intensity. In this case, the only dimensionless parameter toremain constant is the energy (in ω) of field-free ionization (IP /ω) or excitation (En /ω). Additionally, the resonance condition Eq. (3.1) is achieved at different times with a dependence on intensity√ that is not linear (namely, ∝ acos I ). Again, in Fig. 3.7(b) we see good agreement with therising 9ω → 3p resonance delayed as above. The falling resonance appears to have the correct
shape, but appears slightly later than a delay of td = 0.5 cycles.
Next, we scan the driving wavelength (results in Fig. 3.8); the conclusion is similar, though
the width of the scan window is larger, so we see more different features. Generally, radiation
features appear at certain times (0-2 o.c. after the peak of the pulse) then a split occurs atlarger wavelengths. The dependence is clearly nonlinear (and without a familiar analytic form; see
Eq. (3.1)). Usually – but not always – the falling-edge radiation is somewhat more intense than therising-edge. Again, the best/only resonance candidate for the emitted radiation for wavelengths
near 800 nm is 9ω → 3p, with fairly good agreement (incl. offset). For other wavelengths, somefeatures show similar bowl-like structure, while others have a zig-zag pattern of radiation time w.r.t
wavelength.
An additional parameter accessible for variation in computation but, in general, not in experiment is the atomic potential. In Sec. 3.8.1, we will consider helium atoms in the single-activeelectron approximation, but this introduces many differences since the form of the potential changes(affecting different states differently) and equivalent dimensionless parameters require moving to adifferent wavelength and intensity region. Instead, here we consider continuous variations from thehydrogen atom case by altering the nuclear charge, Z, without other modifications. This results inthe single effect that all transition energies scale by Z 2 . We expect 1s →(field-free-)continuum and1s → 3p to have the greatest impact here, but cannot necessarily parse out the consequences of68

Figure 3.5: Temporal profile of radiation within 1ω of the 9th harmonic as a function of pulse
duration in atomic hydrogen for an 800 nm, 3.8 × 1013 W/cm2 pulse. For each duration the results
have been normalized to a maximum of 1 for sake of comparison. A vertical dashed line denotes
the peak of the laser and a horizontal dashed line identifies the pulse parameters used as a base inthis Chapter (i.e., 20-cycle duration).

69

Figure 3.6: Same as Fig. 3.5, but overlaid with lines denoting: (top) resonance times for 7- (magenta) and 8- or 9- (white) photon resonances (there are no resonances with a smaller or larger
number of photons involved for these parameters) for 2s − 6p, and (bottom) 9ω → 3p resonances
delayed by 0.5 (solid), 1.8 (dashed), and 3.0 (dotted) optical cycles.

70

Figure 3.7: Temporal profile of radiation within 1ω of the 9th harmonic as a function of driving
pulse intensity in atomic hydrogen for an 800 nm, 20 o.c. pulse. For each intensity the results
have been normalized to a maximum of 1 for sake of comparison. A vertical dashed line denotes
the peak of the laser and a horizontal dashed line identifies the pulse parameters used as a basein this Chapter (i.e., 3.8 × 1013 W/cm2 intensity). Panels show identical radiation data, but thelower is overlaid with lines denoting 9ω → 3p resonances delayed by 0.5 (solid), 1.8 (dashed), and
3.0 (dotted) optical cycles.

71

Figure 3.8: Temporal profile of radiation within 1ω of the 9th harmonic as a function of driving
wavelength in atomic hydrogen for a 20 o.c. pulse at 3.8 × 1013 W/cm2 . For each wavelength theresults have been normalized to a maximum of 1 for sake of comparison. A vertical dashed line
denotes the peak of the laser and a horizontal dashed line identifies the pulse parameters used as abase in this Chapter (i.e., 800 nm wavelength). Panels show identical radiation data, but the lower
is overlaid with lines denoting 9ω → 3p resonances delayed by 0.5 (solid), 1.8 (dashed), and 3.0(dotted) optical cycles.

72
other state shifts. We note that the field-shifted transition energies do not scale in this same way
since UP does not depend on the potential of the target.
The results of the respective scan are shown in Fig. 3.9, and we see a variety of features quite
similar to those present in the wavelength scan in Fig. 3.8. Again, there are a number of “bowl”
features which match the shape of resonance lines, as well as a few zig-zags that do not. The initial
parameters (Z = 1, λ = 800, and I = 3.8 × 1013 W/cm2 ) all seem to place us just at the transition
point between dominance of two different features, associated with the f alling 9ω → 3p + 0.5 o.c.
and the rising 9ω → 3p + 3.0 o.c. delayed resonance lines.
From these scans, we conclude the identified resonances are correct, and the delays at least
are approximately correct (the falling-side delay consistently appears to be slightly longer than 0.5
cycles, but for now we leave the rising/falling delays equal). Further off of the base parameters,some of the radiation diverges slightly from the 9ω → 3p predictions, while other features clearly
do not correspond to it at all. The deviations are explored in Sec. 3.6, while other resonances may
be identified using the process outlined above for different base scan intensities.

3.4

Radiation at other driving intensities
While Figs. 3.2 and 3.7 show how the generated radiation varies close to 3.8 × 1013 W/cm2 ,it is also valuable to consider a broader range in intensities. In Fig 3.10, we plot the spectrum for
driving intensities ranging from 0.5 − 9.0 × 1013 W/cm2 . From this scan we identify and select fiveintensities of particular interest (vertical white lines):• 2.0 × 1013 W/cm2 : Here are radiation features near 9ω similar to the case of 3.8 × 1013
W/cm2 , studied previously. However, there is additional strong radiation in the region of
the Rydberg state transition energies.
• 3.0 × 1013 W/cm2 : A “normal” case where the harmonic radiation at 9ω clearly dominates
over any off-harmonic radiation.
• 3.8 × 1013 W/cm2 : The case we have investigated most closely, where radiation just below

73

Figure 3.9: Temporal profile of radiation within 1ω of the 9th harmonic as a function of atomic
nuclear charge Z in a hydrogen-like atom for an 800 nm, 20 o.c. pulse at 3.8 × 1013 W/cm2 . For
each Z the results have been normalized to a maximum of 1 for sake of comparison. A verticaldashed line denotes the peak of the laser and a horizontal dashed line identifies the pulse parameters
used as a base in this Chapter (hydrogen, i.e., Z = 1). Panels show identical radiation data, but
the lower is overlaid with lines denoting 9ω → 3p resonances delayed by 0.5 (solid), 1.8 (dashed),and 3.0 (dotted) optical cycles.

Nγ ω → |ni
9ω → 3p
9ω → 3p
9ω → 3p
10ω → 4s
10ω → 3s
9ω → 2p
7ω → 2p
9ω → 6p(or 4d)(or 3d)(or 6f or 6h)

td (o.c.)
0.5
1.8
3.0
0.5
1.2
0.8
0.8
0.8

line style
white
white
white
magentaredcyanyellow
green

74
solid
dasheddottedTable 3.1: Identified resonances, delays, and line styles. We note that due to l-independent energies
in hydrogen atoms, some of the resonances can be ascribed to different states with the same l parity,
namely: 3s and 3d, 4s and 4d, 6p and 6f or 6h.

and above 9ω dominates over that at the harmonic.
• 5.6 × 1013 W/cm2 : Similar features to the radiation generated at 3.8 × 1013 W/cm2 , butwith strong and clearly-separated sidebands well above and below 9ω.
• 7.0×1013 W/cm2 : Here, again the radiation emitted at 9ω dominates, but weaker sidebands
are still present.
For each of these base intensities we have performed the same process to identify resonancesand delays. As above, we generally have considered – in order – scans across pulse duration,
driving intensity, and driving wavelength. From each scan, we can eliminate potential resonances
which do not match the shape of the radiation curves, while estimating delays for the remaining
possibilities. Through this, we identified resonances and delays, accounting for most of the clear
radiation structures. In the Figures below, we plot all identified resonance-delays, labeled byresonance and denoted by color defined in Table 3.1. As mentioned above, 9ω → 3p is the only
resonance for which we have identified multiple delays, differentiated by line style.
First, in Figs. 3.11 and 3.12, we show scans of pulse duration for each chosen intensity. As
above, all resonance conditions are straight lines versus pulse duration, so it is simple to identify
radiation connected to resonances. As in Fig. 3.5, which radiation feature is dominant often depends
on pulse duration. However, we are not able to identify a systematic rule or trend for this variation.
We note a plotting artifact visible in some figures (esp. in the lower part of Fig. 3.11’s top panel):

75

Figure 3.10: HHG spectra as a function of peak intensity for radiation during the pulse. Resultsfor each intensity have been normalized to the maximum signal in the region between 8.75ω and9.5ω. Red dashed lines denote field-free excited state transition energies.

76
As a function of duration, it appears that radiation rapidly transitions between peaks and nodes.
This is a result of relatively coarse sampling, not a physical effect – radiation is slice-normalized foreach sampled duration, but smoothing in the Figure interpolates lower values for pixels between
actual data.
For 2.0 × 1013 W/cm2 (Fig. 3.11 top), a duration-independent radiation feature just after
the peak of the pulse dominates. This feature corresponds to the “normal” harmonic emission at
9ω, without resonant enhancement. However, for all durations, there are two additional durationdependent radiation features visible, which we identify to correspond to rising and falling resonances. Specifically, the 7ω → 2p + 0.8 o.c. (yellow) and 9ω → 6p + 0.8 o.c. (green) delayedresonances both agree well with these features on both the rising and falling sides (incl. across otherscans, see below). From the present results, we can however not determine conclusively which one(or both) is the origin of the generated radiation. We note that the same resonances occur at all
higher intensities, but are in the rising and falling edge of the pulse far from the peak and in thosecases they do not appear to provide significant radiation.
In the spectra at 3.0 × 1013 W/cm2 (Fig. 3.11 bottom) the dominant radiation occurs justafter the peak of the pulse for durations up to about 15 o.c., then it shifts to earlier times withincrease of the pulse duration. This suggests that resonances do not play a strong role for shorter
pulses at this intensity. We have not definitively identified a resonance associated with the longerpulse radiation. There is some agreement (across different scans) with a much larger delay of
the resonances identified at the lower intensities (7ω → 2p and 9ω → 6p), specifically 4.3 o.c..
Unfortunately, the agreement is generally good only in small portions of each scan near the baseparameters (800 nm, 20 o.c., 3.0 × 1013 W/cm2 ). Unlike the delayed resonances listed in Table 3.1,
here the features do not correspond with radiation in different and broader scans (for example, inFig. 3.12, we do not identify radiation features with the same slopes as the green/yellow lines atlater delays).
For 5.6×1013 W/cm2 (Fig. 3.12 top), there are two primary radiation structures corresponding
to the 10ω → 4s + 0.5 o.c. delayed resonance (magenta), but other structures are also present.

77Firstly, unlike for 3.8 × 1013 W/cm2 , there are particular times (separated by 0.5 o.c.) whenradiation is emitted regardless of pulse duration. The intensity of this radiation increases strongly
when it coincides with resonance lines (esp. 10ω → 4s). This suggests a combined influence of the
field and envelope, which we will further investigate in Sec. 3.5. We note also the presence of otherresonance lines. Rising 9ω → 3p + 1.8 o.c. may be responsible for broadening and “smearing” therising radiation feature, and weaker features line up with other resonances.
Finally, at 7.0 × 1013 W/cm2 (Fig. 3.12 bottom) we can see a single strong radiation featurecorresponding to rising 10ω → 3s + 1.2 o.c. (red) and 9ω → 2p + 0.8 o.c. (cyan) resonances, but agreat deal of weaker radiation features are also present. Similarly to 5.6 × 1013 W/cm2 , radiationis enhanced at particular times, but the effect is not as sharp.
In Fig. 3.13, we present the temporal profile of radiation for the full intensity range. Before
considering individual radiation features, we may make some general observations. As for intensities
near 3.8 × 1013 W/cm2 , the radiation appears in several “bowl” shapes centered just after the peakof the pulse. For lower intensities ( / 5 × 1013 ), the structure of radiation tends to be smootherand be strongest within the first three cycles after the peak of the pulse. At the higher intensities,
other structures are visible (e.g., discrete radiation times independent of intensity). Similarly, weobserve more complexity in radiation structures occurring after the peak of the pulse compared tothose before the peak.
With good agreement we have included identified resonances and delays in Fig. 3.13 (bottom
panel), which clearly explain the majority of radiation features. In regions where two or more
delayed resonances overlap, the radiation features broaden, but the number of resonances does notappear to predict which radiation will dominate (cf. the 5 − 6 × 1013 W/cm2 range).
Scanning the wavelength for each of the four intensities (Figs. 3.14 and 3.15) serves as animportant additional check in resonance identification since the timing has a more complex dependence than in either the intensity or duration scans. In general, the agreement is quite good withina range of 50 nm centered about 800 nm. Beyond this range, there are radiation features whichappear to correspond to unidentified resonances. In other cases identified resonances seem to have

78

Figure 3.11: Temporal profile of radiation within 1ω of the 9th harmonic as a function of pulse
duration in atomic hydrogen for an 800 nm pulse at (top) 2.0×1013 W/cm2 and (bottom) 3.0×1013
W/cm2 . For each duration the results have been normalized to a maximum of 1 for sake ofcomparison. A vertical dashed line denotes the peak of the laser and a horizontal dashed line
identifies the pulse duration used as a base in this Chapter (i.e., 20 o.c.). Lines are drawn for
delayed resonances as outlined in Table 3.1.

79

Figure 3.12: Same as Fig. 3.11, but for (top) 5.6 × 1013 W/cm2 and (bottom) 7.0 × 1013 W/cm2 .

80

Figure 3.13: Temporal profile of radiation within 1ω of the 9th harmonic as a function of driving
pulse intensity in atomic hydrogen for an 800 nm, 20 o.c. pulse. For each intensity the results havebeen normalized to a maximum of 1 for sake of comparison. A vertical dashed line denotes thepeak of the laser and horizontal dashed lines identify the intensities selected above. (lower plot)
Lines are drawn for delayed resonances as outlined in Table 3.1.

81
the correct behavior with changing wavelength, but occur at larger wavelengths than the corresponding radiation (i.e., lines appear at vertical offsets from radiation peaks). Identifying these
additional resonances would require more simulations (e.g., scanning various parameters around abase pulse of 700 nm and 2.0 × 1013 W/cm2 or 875 nm and 5.6 × 1013 W/cm2 ). Deviations in thealready-predicted resonances are more interesting, potentially highlighting approximations in ourmodel; we will therefore investigate these in Sec. 3.6.

3.5

CEP and the trajectory modelPreviously, Camp et al. [100] ascribed similar features to resonantly-initiated trajectories,

and described their free parameters as falling within the range of predictions from the group’s
earlier semi-classical trajectory calculations in [41]. Using a classical trajectory code based off thedescription in this paper, we however did not find consistent initial conditions leading to trajectoriesof the correct durations (0.5, 1.8, 3.0 o.c.) for our radiation at 3.8 × 1013 W/cm2 .
We have also tested the semi-classical model based on the following considerations. If semiclassical trajectories are involved in the observed radiation, we would expect a clear dependenceof the radiation features on the carrier-envelope phase (CEP) since in any semi-classical trajectorymodel, field values between ionization and recombination times will determine the travel time andrecombination energies possible. In “normal” HHG, the timing of radiation at a given energy
is constrained by the phase of the laser at time of ionization. Further, in the tunneling regime,
the field strength at moment of ionization determines both the ionization probability and theinitial parameters of trajectories (location of the tunnel exit). In a theoretical study where timeof ionization was gated by application of a second, VUV laser, field values at the moment ofionization continued to strongly impact the spectral and temporal character of HHG [102]. Resultsfrom a generalized semi-classical trajectory model mimicking trajectories initiated by multiphoton
ionization showed a dependence at least on the direction of the field, differentiating “uphill” and“downhill” trajectories [41].
In Fig. 3.16, we investigate the CEP dependence for two intensities – (top panel) 3.8 ×

82

Figure 3.14: Temporal profile of radiation within 1ω of the 9th harmonic as a function of driving
pulse wavelength in atomic hydrogen for a 20 o.c. pulse at (top) 2.0 × 1013 W/cm2 and (bottom)
3.0 × 1013 W/cm2 . For each wavelength the results have been normalized to a maximum of 1 forsake of comparison. A vertical dashed line denotes the peak of the laser and a horizontal dashedline identifies the laser wavelength used as a base in this Chapter (i.e., 800 nm). Lines are drawn
for delayed resonances as outlined in Table 3.1.

83

Figure 3.15: Same as Fig. 3.14, but for (top) 5.6 × 1013 W/cm2 and (bottom) 7.0 × 1013 W/cm2 .

84
1013 W/cm2 , and (bottom) 5.6 × 1013 W/cm2 – and discover an interesting contrast. For the
lower intensity, radiation appears to be completely independent of CEP, consistently agreeing withdelayed resonance times. This suggests that the semi-classical trajectory model is not applicablehere. However, for the higher intensity, we see a clear effect, where the timing of the dominantradiation peak shifts linearly with CEP. As in earlier scans, the radiation is strongest at discrete
times (every half-cycle) especially when close to the timing of the delayed resonance (f alling 10ω →4s + 0.5 o.c.). While in previous scans, these discrete times were fixed, here they shift in the same
way as the electric field, with a slope of 0.5 cycles per π-phase shift. The consistent timing ofradiation – independent of resonance enhancement to its intensity – suggests similarity between
multiphoton- and tunnelling-initiated trajectories at this intensity. Further study at the higher
intensities will be necessary to verify this observation.
Thus, at lower intensities (3.8 × 1013 W/cm2 , as well as 2.0 and 3.0 × 1013 W/cm2 , notshown), this resonant radiation appears to be fully controlled by the laser envelope independent of
the specific phase of the electric field. In contrast, at higher intensities (5.6 and 7.0 × 1013 W/cm2 ),there is an interplay between the envelope effect (resonance time) and the field effect (trajectories
and/or time of tunnel ionization). We however also note that the Keldysh parameter, typicallyused to separate multiphoton (γ > 1) from tunneling (γ < 1) behavior, falls in the multiphoton
regime for all intensities – γ = 1.7 at 3.8 × 1013 W/cm2 and γ = 1.4 at 5.6 × 1013 W/cm2 .
In Fig. 3.17, we apply the frequency-offset prediction of [100] (Eq. (3.2)) with α = 2π forscans of intensity and wavelength. We observe good agreement for certain features (e.g., red-shifted
radiation near 6 × 1013 W/cm2 ). However, the model does not well account for the majority of
blue-shifted radiation – especially that near 2 × 1013 W/cm2 – which displays a clear left-opening
shape in contrast to the right-opening predictions.

3.5.1

Varying pulse envelope

In seeking to parse the effect of the field at the time of resonance versus the time of radiation,
we additionally vary the envelope of the electric field. Specifically, we steepen either the rising or

85

Figure 3.16: Temporal profile of radiation within 1ω of the 9th harmonic as a function of driving
pulse CEP in atomic hydrogen for an 800 nm, 20 o.c. pulse at (top) 3.8×1013 W/cm2 , and (bottom)
5.6 × 1013 W/cm2 . A vertical dashed line denotes the peak of the laser and horizontal dashed lines
identify the CEP shown in other scans (φ = π/2). Lines are drawn for delayed resonances as
outlined in Table 3.1.

86

Figure 3.17: HHG spectra as a function of (top) peak intensity at λ = 800 nm and (bottom)
driving wavelength at I = 3.8 × 1013 W/cm2 , in atomic hydrogen for a 20 o.c. pulse. We applieda Blackman window to dipole acceleration, and results for each intensity (top) or each wavelength(bottom) have been normalized to the maximum signal in the region between 8.75ω and 9.5ω.
Overlaid curves denote the ∆ω prediction, Eq. (3.2), for (top) delayed resonances as outlined inTable 3.1 and (bottom) 9ω → 3p resonances delayed by 0.5 (solid black), 1.8 (dashed black), and
3.0 (dotted black) optical cycles.

87
falling half of the envelope by using higher powers of the sine function – sin4 , sin6 or sin8 rather
than the standard sin2 . We note that these even-integer powers retain necessary continuity anddifferentiability. Varying the envelope serves also as a tool for identifying resonance delays, andwas valuable in the initial identification of the 9ω → 3p resonance and delays.
In Fig. 3.18, we plot temporal profiles for various envelopes. We note that plotting of thisdiscrete scan differs from other Figures, with results for each envelope appearing as an independent
horizontal bar. We first observe that the most intense (last) feature is minimally affected (timing
or intensity) by modification of the rising side of the pulse, but shifts earlier when the falling sidesteepens. Likewise, the earliest feature is unaffected by the falling side (as required by causality)
but shifts later as the rising side steepens. Considering the second feature – which occurs right at
the peak of the pulse for sin2 ramp-up – we see it shifts in time similarly to the first, with a slightincrease in intensity for envelopes with steeper rise. Finally, the third feature’s (most noticeable forthe original sin2 envelope) position does not change when the back of the pulse steepens, but shifts
later (merging with the fourth) as the front of the pulse steepens. The intensity of the third featuredoes not appear to change greatly, but it is difficult to distinguish from the fourth feature for bothsteepest envelopes. In summary, the character of the earliest three radiation features appears to becontrolled by the rising side of the pulse, including the feature occurring a full cycle after the peak.

3.6

Deviations from predicted resonance times
In various scans above, we noted that the observed radiation features generally fit well withpredictions of delayed resonances, but that the two sometimes disagreed somewhat for parts of thescan. Thus, we sought to further test the consistency of these delays. A simple way to do this is tochoose a particular resonance and fix the resonance time, Tr (in cycles from peak) across scans. One
would then expect that the respective radiation feature will occur along vertical lines on temporalprofile scans, thus clearly highlighting potential deviations and trends.

88

Figure 3.18: Temporal profile of radiation within 1ω of the 9th harmonic for a selection of electricfield envelopes (horizontal bars) in atomic hydrogen for an 800 nm, 20 o.c. pulse at 3.8 × 1013
W/cm2 . We note that the results have not been normalized. A vertical dashed line denotes thepeak of the laser. The center bar corresponds to the sine-squared envelope used elsewhere in thisthesis. In the lower panel, lines are drawn for delayed resonances as outlined in Table 3.1.

89
To this end, we again consider the resonance condition


I
Tr π
4
cos ±
= Nγ − Z 2 (En − E0 ),
4ω 2
Nτ

(3.3)

as in Eq. (3.1). Here, we however fix Ei as the hydrogen atom energies but allow variation of the
Coulomb charge, Z. In order to keep Tr fixed, it is required to vary multiple parameters together.
Still, many different variations are possible within the parameter space. To constrain our investigation, we vary pairs of continuous parameters: pulse duration in cycles (Nτ ), Coulomb charge (Z),wavelength (λ ∝ 1/ω), and intensity (I). We did not vary the envelope shape since only discrete
sin2n envelopes are viable, and excluded the CEP as a parameter since it does not affect resonancetime or – as shown previously – radiation at low peak intensities. These considerations result in six
distinct scans for each resonance (defined by fixing Nγ and En ), where we scan one parameter overa certain range while changing the second parameter accordingly to keep Tr constant, as describedin Table 3.2. The parameter range is chosen such that the HHG cut-off at IP + 3.17UP remains
consistently larger than 10ω.

3.6.1

Deviations in radiation features around 9-photon resonance to 3p state

To investigate the 9ω → 3p resonance, we consider scans in the parameter subspace definedby Nγ = 9, En = E|3pi = −0.056 a.u., and Tr = 2.0 o.c., which includes our primary case ofinvestigation analyzed in detail above, i.e., I0 = 3.8 × 1013 W/cm2 , λ0 = 800 nm (ω0 = 0.057 a.u.),Nτ,0 = 20 o.c., Z0 = 1. As mentioned above, by scanning along a variety of curves with constant Tr ,we intend to identify the stability of radiation features as well as the potential effects of the three
dimensionless parameters. In particular, all scans shown in Table 3.2 have the common signaturethat z = UωP decreases or remains fixed along the scan, while there are different trends for IωP andUP
IP .

First, we demonstrate the importance of dimensionless parameters by considering a threeparameter scan which keeps simultaneously constant all of:

UP I P UP
ω , ω , IP , and Tr

for each of theresonances. We do this by scanning Z (from 0.9 to 1.1) and fixing λ = 800
nm (ω ∝ Z 2 ) and
Z2

1
1

10, 10.5...30
10, 10.5...30

10, 10.5...30

20Nτ − ZNτ − IBase

1

En −E0INγ ω0 − 02 cos4
4ω0

1

τ



± TNr π4ω 2 (Nγ ω−Z02 (En −E0 ))Tr π
cos4 ± N

750, 755...850800

3.8 × 1013

1.47

2.8 & 1.3

τ



Tr π
cos4 ± N
τ,0I0 4  Tr π 
cos ± N800

1.7 & 1.4

3.8 × 1013

numericsolve

1.5

3.8 × 1013

800

τ,0

2.1 & 0.9

2.4 & 0.5τ,04ω02 (Nγ ω0 −Z 2 (En −E0 ))Tr π
cos4 ± N

800

z = UωP
2.2 & 1.0I
W/cm2
3.8 × 1013

λnm
numeric
solve

8.78

8.8

9.3 & 8.7

9.4 & 8.6

8.2 % 9.3

7.9 % 9.7

8.1 % 9.2

IP
ω

0.167

0.32 & 0.15

0.19 & 0.16

0.16 % 0.17

0.25 & 0.09

0.30 & 0.05

0.27 & 0.10

UP
IP

10.24

11.6 & 10.1

11.0 & 10.1

10.9 & 10.1

10.3 & 10.2

10.3 & 10.2

10.3 & 10.2

IP +UPω

Table 3.2: Classification of two-parameter scans keeping 9ω → 3p resonance time Tr fixed. For each scan, two parameters are fixed attheir base values and two are varied: one parameter is given as a range of equally-spaced values, while the other is determined by solving
Eq. (3.3) analytically or numerically. The final four columns describe the range and trends of dimensionless parameters across the scans.

Nτ − λ

20

λ−I
r

20

Z −I

0.90, 0.905...1.100.95, 0.9525...1.05

20

Z −λ

Z

Nτ
o.c.

Title

90

91
I = Z 6 ∗ 3.8 × 1013 W/cm2 . With this variation for the TDSE in length gauge:


∂
∇2
1 ~
i ψ(~r, t) = −+ − + E(t) · ~r ψ(~r, t),
∂t
2
r

(3.4)

√
with E → E0 I = E0 Z 3 and − 1r → − Zr , we get


∇2
Z
∂
3 ~
+ − + Z E0 (t) · ~r ψ(~r, t),i ψ(~r, t) = −
∂t
2
r
0

(3.5)

0

which is equivalent to Eq. (3.4) under the transform t → Zt 2 and ~r → ~rZ :
i

∂
ψ
∂t

 0  0  0 0 
 0 0  
t
∇2
~r
~r t
~r t
Z2
·
ψ
.
, 2 = −
+ − 0 + Z 3 E~0
,
Z Z
2
r
Z2
Z
Z Z2

(3.6)

Thus, the system has merely been rescaled (please note that the derivation is much the same in thevelocity gauge). Consistently, the results in Fig. 3.19 demonstrate that a simultaneous variationof wavelength across 660 − 988 nm and intensities from 2.0 − 6.7 × 1013 W/cm2 along with avariation of the Coulomb charge Z produces the same radiation (up to a factor). Beyond puremathematical interest, this highlights a utility of dimensionless parameters over traditional ones.
This also happens to be a fair test of the numerical stability of our code – we did not alter our gridor other numerical parameters and yet achieved consistent observables (incl. those not shown suchas full complex spectrum through the cut-off, total ionization, and excited state population after
the pulse).
Now, we turn to scans of more practical interest for the present purpose. We begin byconsidering the first three types of scans in Table 3.2. In each of these scans, two of the threeparameters Z, λ, and I are varied. This results in a change of all three dimensionless parameters,
however the trends for the dimensionless parameters are the same across the scans ( UωP and UIPP
decreasing, IωP increasing). Since the trends are opposite for IωP (increasing) and UωP (decreasing),PP
the energy (in ω) required for field-shifted ionization, IP +U, and for resonant excitation, En +U,ωω

at the peak of the pulse remains almost the same, up to about 0.1.
The results of the respective scans are shown in Fig. 3.20 and surprisingly show similarradiation trends in each of these scans. To avoid implying a specific association between trends anda particular parameter (e.g., Z or UωP ), we discuss the results in terms of progress along a scan or

92

Figure 3.19: Temporal profile of radiation within 1ω of the 9th harmonic as a function of atomic
nuclear charge Z in a hydrogen-like atom while also varying wavelength (λ = 800nm) and intensity
Z2(I = Z 6 ∗ 3.8 × 1013 W/cm2 ) to keep resonance times fixed. A 20 o.c. pulse was used throughout.
For each Z the results have been normalized
to a maximum of 1 for sake of comparison (note:√
spectral intensity here scales roughly ∝ Z). A vertical dashed line denotes the peak of the laser
and a horizontal dashed line identifies the pulse parameters used as a base in this Chapter (i.e.,
Z = 1, I = 3.8 × 1013 W/cm2 , λ = 800 nm). Overlaid with lines denoting resonance times for 7(magenta) and 8- or 9- (white) photon resonances.

93
location in the associated figure, where “up” is further up on the respective plot, and corresponds
to an increase of the first scanned variable. We first observe a transition between dominance of theouter radiation features (earlier associated with rising and f alling 9ω → 3p + 0.5 o.c.) and theinner ones (rising 9ω → 3p + 1.8 o.c. and rising 9ω → 3p + 3.0 o.c.) near the base parameters,
which are denoted by the horizontal white line in each plot. Above this point (e.g., Z > 1 for the
top panel), the outer features dominate, while the inner features dominate below it. Additionally,
the outer radiation features shift inwards significantly toward the top of the scans, tending to occur
closer to the peak of the pulse. At the same time these features decrease in magnitude, as canbe best seen in Fig. 3.21, where we show the same data as in Fig. 3.20 (top panel) but this time
un-normalized. Near the bottom of the scans, the inward trend is smaller, but still noticeable. The
inner radiation also decreases in magnitude, even more rapidly, when moving up the scan, but thetiming is more consistent. The peak of the second feature (rising 9ω → 3p + 1.8 o.c.) appears
slightly earlier in time as one moves down in the scan. Actually, the prediction for the resonance time
is best matched slightly below the base parameters. However, the un-normalized plot (cf. Fig. 3.21)
suggests that the radiation feature is indeed more widening in time than shifting. The third feature
stays well-matched with the prediction for all portions of the scans, as long as it is distinguishable.
Other identified resonance delays do not appear to play a role, nor do resonances not previously
associated match the shape of radiation curves. From these results, we can again definitively statethat radiation behavior depends not on traditional parameters (I, λ, Z) but rather on an alternate
set such as the dimensionless parameters of Table 3.2. Despite very different traditional parameterchanges across the three scans, radiation looks much the same for each scan. Unfortunately, we
cannot tease out which dimensionless parameter(s) predict behavior – UωP , UIPP , and IωP have similar
trends for all three scans. Further, the Z − I scan has the widest range for each parameter,
consistent with the impression of Fig. 3.20). Thus, the loss of agreement with resonance predictions
may be attributed to: perturbative behavior (decreasing UωP ), greater multiphoton dominance overtunneling (decreasing UIPP ), more photons involved in field-free interactions (increasing IωP ), or somecombination.

94

Figure 3.20: Temporal profile of radiation within 1ω of the 9th harmonic for scans of: (top) atomic
nuclear charge Z and wavelength λ, (middle) Z and intensity I, and (bottom) λ and I. For allscans, the 9ω → 3p resonance time is kept fixed. A 20 o.c. pulse was used throughout. For eachelement of the scan (i.e., horizontal slice), the results have been normalized to a maximum of 1 forsake of comparison. A vertical dashed line denotes the peak of the laser and a horizontal dashedline identifies the pulse parameters used as a base in this Chapter (i.e., Z = 1, I = 3.8 × 1013
W/cm2 , λ = 800 nm). Lines are drawn for delayed resonances as outlined in Table 3.1.

95

Figure 3.21: Same as top panel of Fig. 3.20, but with no normalization for each Z and on alogarithmic scale. Unscaled versions for other panels show similar trends, except with an overall≈ I 4 scaling of the radiation magnitude.

96
For the remainder of the scans, the primary scanned variable is pulse duration in optical
cycles, with I, λ, or Z adjusted to keep Tr = 2. For these scans the set of three dimensionlessPPparameters vary differently (see Table 3.2), though IP +U
and En +U
do have similar trends.
ωωAdditionally, it is worth considering Nτ itself as a meaningful dimensionless parameter describing
the (inverse) spectral bandwidth and the inter-cycle envelope stability.
In Fig. 3.22, there is a similar transition in dominance of the features across the base parameters that we have noticed in the first three scans. The outer (inner) radiation dominates at
the top (bottom) of the scans. In this case, the absolute intensity of the outer features increases
up the scan, while the inner ones decrease (cf. un-normalized plot Fig. 3.23). The peaks of theouter features (rising and f alling 9ω → 3p + 0.5 o.c.) show a clear ‘V’-like trend away from thepeak of the pulse. This follows the shape of the lower-intensity resonance lines such as 9ω → 4p,
possibly suggesting that the 3p state shifts slightly less than UP and the scans are highlighting thisdeviation. On the other hand, the results in Fig. 3.23 de-emphasize the shifting of the peaks andagain more clearly show that the outer features widen away from the prediction. The timing of
the second radiation feature (rising 9ω → 3p + 1.8 o.c.) appears consistent with the prediction– constant. The third feature (rising 9ω → 3p + 3.0 o.c.) – which dominates the radiation for
shorter-than-base pulses – is more complicated, with the peak appearing to oscillate slightly earlier
and later (overall slightly earlier than prediction). At the bottom of the scans, radiation diverges
from the 9ω → 3p resonances and appears to correspond to other resonances (namely, 10ω → 4s,
10ω → 3s and/or 9ω → 2p). For the three scans, we contrast the similarity of radiation structures
to the differences in trends of the dimensionless parameters UωP , UIPP , and IωP . Therefore, it seems
likely that the effects are a result of increasing pulse duration itself (potentially via extending duraP
tions of resonance) and altering the resonances occurring during the pulse (for which IP +U
servesω

as a proxy).
Continuing the discussion of the first three scans, we note that values of IωP ≈ 9.3 foundat the top of Fig. 3.20 do not appear to disrupt resonant behavior near the bottom of the upper
panels in Fig. 3.22. However, low values of UωP (/ 1) and UIPP (/ 0.1) are not included in these latter

97
scans, so we may make no further conclusions to their roles. Overall, the two trios of scans seemto explore somewhat orthogonal spaces.
When varying more than two parameters simultaneously, it is possible to construct scanswith other trends of dimensionless parameters (e.g., keeping UIPP constant while varying UωP andIPω ), and we recommend this for future investigation.

3.7

Short and long pulses
In prior results, we have focused on pulses with durations near 20 cycles. It is valuable now

to briefly discuss effects when the driving pulse has either much shorter or much longer duration.
For short pulses, the envelope variation may be significant even within a single cycle, and the“average” energy UP (t) becomes less sensible. On the other hand for long pulses, the relative
intensity stability means both that resonances occur over longer stretches of the pulse and that theenvelope does not vary significantly over the period between resonance and radiation.
For a 4-cycle pulse of 3.8 × 1013 W/cm2 , we perform a scan of CEP (Fig. 3.24). We observe
radiation at CEP-dependent times, suggesting trajectory-like behavior as seen only at the higher
intensity in Sec. 3.16. However, this radiation shows no enhancement when it coincides withpredicted radiation times. Instead, radiation is strongest just at and after the peak of the pulse, as
expected for “normal HHG.”In Fig. 3.25, we extend the scan of Fig. 3.6 through 50 o.c.. For pulses longer than about 30
o.c., the radiation predictions begin to break down. First, the majority of radiation occurs within
the first two cycles after the laser peak, not matching any of the previously-identified resonancedelays. Secondly, the earliest radiation feature begins to shift earlier than the rising 9ω → 3p +0.5 o.c. delayed resonance. Agreement with other resonances would require large delays (greaterthan 3 cycles) and initial investigations find poor agreement.

98

Figure 3.22: Temporal profile of radiation within 1ω of the 9th harmonic for scans of pulse duration,Nτ , in optical cycles and (top) atomic nuclear charge Z, (middle) wavelength λ, and (bottom)
intensity I. For all scans, the 9ω → 3p resonance time is kept fixed. For each element of the scan(i.e., horizontal slice), the results have been normalized to a maximum of 1 for sake of comparison.
A vertical dashed line denotes the peak of the laser and a horizontal dashed line identifies the pulse
parameters used as a base in this Chapter (i.e., Z = 1, Nτ = 20 o.c., I = 3.8 × 1013 W/cm2 , λ = 800
nm). Lines are drawn for delayed resonances as outlined in Table 3.1.

99

Figure 3.23: Same as middle panel of Fig. 3.22, but with no normalization for each Nτ and on alogarithmic scale. Unscaled versions for other panels show similar trends, except with an overall≈ I 4 scaling of the radiation magnitude.

Figure 3.24: Temporal profile of radiation within 1ω of the 9th harmonic as a function of driving
pulse CEP in atomic hydrogen for an 800 nm, 4 o.c. pulse at 3.8 × 1013 W/cm2 . A vertical dashedline denotes the peak of the laser. Lines are drawn for delayed resonances as outlined in Table 3.1.

100

Figure 3.25: Same as Fig. 3.5, but extending to durations up to 50 o.c..

101

3.8

Appearance in other parameter regimes
While the appearance of this resonance-enhanced behavior requires a great deal of further

study within the relatively narrow parameter space described above, we would emphasize that theeffects are not limited to hydrogen or near-IR lasers. To that end, we present preliminary results
for atomic hydrogen driven by a mid-infrared (1600 nm) laser and helium atoms driven by a laseron the edge close to the UV regime (400 nm). We also recall that Camp et al.’s investigation [100]was for argon atoms driven by a near-IR laser.

3.8.1

400 nm helium

For helium, we perform simulations using the SAE potential presented in Ch. 2, and selecta laser with similar dimensionless parameters to that of the 3.8 × 1013 W/cm2 hydrogen results
above. Namely, we select a wavelength of 400 nm such that UP /ω = 8.1 (vs. 8.8 above), and an
intensity 3.3×1014 W/cm2 such that energy in units of ω remains similar for the HHG cut-off (13.1ω
vs. 13.4ω) and the ponderomotive energy (1.6ω vs. 1.5ω), as well as a similar Keldysh parameter(γ = 1.6 vs. 1.7). We retain the same number of cycles, i.e., 20. The most notable features of theresultant spectrum (Fig. 3.26) are a peak below 7ω (broader and weaker than for hydrogen) anda 9ω radiation suppressed relative to both nearby frequencies and neighboring harmonics, just likein hydrogen. Varying driving laser intensity around 3.3 × 1014 W/cm2 and wavelength around 400
nm, we again see (Fig. 3.27) a clear region where radiation below and above 9ω dominates thatat the harmonic frequency. In Fig. 3.28, we show temporal profiles for scans of (top) duration,(middle) peak intensity, and (bottom) driving wavelength. The structures appear consistent withthose found for hydrogen atoms (Figs. 3.5-3.8).

3.8.2

1600 nm hydrogen

HHG driven by mid-infrared pulses has received growing interest in recent years [103–105].
Due to the scaling of the HHG cut-off (IP +3.17 4ωI 2 ) and lower probabilities of ionization (i.e., higher

102

Figure 3.26: High harmonic spectrum generated in helium, exhibiting a rich structure in the nearthreshold region of the spectrum between harmonics 5 and 12. The helium ionization threshold of
0.94 a.u. = 8.1ω, where ω = 0.114 a.u. corresponds to the central driver wavelength of 400 nm.
Other laser parameters: sin2 pulse of 20 optical cycles with peak intensity of 3.3 × 1014 W/cm2 .
Note that the intensity of the spectrum is scaled by the signal at the fundamental frequency. (cf.
hydrogen results in Fig. 3.1).

103

Figure 3.27: HHG spectra as (top) a function of driving wavelength at I = 3.3 × 1014 W/cm2 and(bottom) a function of peak intensity at λ = 400 nm. Nτ = 20 o.c. for both panels. Results foreach wavelength (top) or each intensity (bottom) have been normalized to the maximum signal
in the region between 8.75ω and 9.5ω. Red dashed lines denote field-free excited state transitionenergies. (cf. hydrogen results in Fig. 3.2 and Fig. 3.3).

104

Figure 3.28: Temporal profile of radiation within 1ω of the 9th harmonic as a function of (top)pulse duration, (middle) peak intensity, and (bottom) driving wavelength. All in helium atom fora pulse of (unless specified) 20 o.c., 400 nm, and 3.3 × 1014 W/cm2 . For each element of the scans(i.e., horizontal slice), the result has been normalized to a maximum of 1. A vertical dashed line
denotes the peak of the laser and a horizontal dashed line identifies the base pulse parameters.

105
saturation intensities), longer driving wavelengths allow for production of bright high-harmonicsupercontinua extending into keV energies [31]. Much of the research on mid-IR HHG focuses onhigh driving laser intensities and the resulting ultrahigh harmonics. Here, we seek to investigate
the near-threshold response, where low photon energies and large UP mean that many states crossthrough resonance.
In the case of a 1600 nm laser, the field-free 1s → 3p transition energy is nearly resonant withthe 13th harmonic (E1s→3p = 13.16ω), and radiation at this energy is consistently more intensethan other harmonics (excl. ω and 3ω). In Fig. 3.29, we show full (top panel) and near-threshold(bottom panel) HHG spectra for two driving intensities. Above IP , the majority of radiation peaks
do not clearly correspond to odd harmonics. Below Ip , harmonic radiation is dominant, but otherstructures are also present, similar to those found with 800 nm pulses.
Fig. 3.30 demonstrates the complex temporal structure of the 13th and 17th harmonics
as a function of peak intensity. We note that many field-shifted resonances occur in this range –overlaying resonance conditions would cover nearly half of the Figure pixels within 3 o.c. of the peak.
For most intensities (above about 1013 W/cm2 ), dominant radiation occurs within about two cyclesof the peak of the pulse. The specific timing, however, is highly dependent on driving intensity,with few consistent trends. Radiation does seem to occur most often centered at 0.5, 1.0, 1.5, and
2.0 cycles – as would be expected for trajectory-driven HHG (see Sec. 3.5) – though we do not
observe patterns as to which of these times will show the strongest radiation.
Considering the lowest intensities, 13th harmonic radiation shows cleaner structure, andwe investigate this in Fig. 3.31. Here, the majority of radiation occurs before the peak of thepulse, following the shape of various rising-side resonances (esp. 17ω → 5p (or 5f ) and 16ω →3s (or 3d)). Additional weaker (and later) structures show similar trends, suggesting the role ofmultiple resonances and/or delays.

106

Figure 3.29: High harmonic spectra generated in atomic hydrogen. (top) Broad plateau extending
to cut-offs (near 27 and 67ω, denoted by dashed lines) with many peaks appearing offset fromN ω; (bottom) near-threshold region exhibiting a rich structure. The hydrogen ionization thresholdof 0.5 a.u. = 17.5ω, where ω = 0.028 a.u. corresponds to the central driver wavelength of 1600
nm. Other laser parameters: sin2 pulse of 10 optical cycles with peak intensities of 0.75 × 1013
and 4.95 × 1013 W/cm2 . Note that the intensity of the spectra are scaled by the signal at thefundamental frequency.

107

Figure 3.30: Temporal profile of radiation within 1ω of the (top) 13th and (bottom) 17th harmonics
as a function of driving pulse intensity in atomic hydrogen for a 1600 nm, 10 o.c. pulse. For eachintensity the results have been normalized to a maximum of 1 for sake of comparison. A verticaldashed line denotes the peak of the laser and horizontal dashed lines identify the intensities inFig. 3.29.

108

Figure 3.31: Temporal profile of radiation within 1ω of the 13th harmonic as a function of driving
pulse intensity in atomic hydrogen for a 1600 nm, 10 o.c. pulse. For each intensity the results
have been normalized to a maximum of 1 for sake of comparison. A vertical dashed line denotes
the peak of the laser and a horizontal dashed line identifies the lower intensity in Fig. 3.29 (i.e.,
0.75 × 1013 W/cm2 ). Panels show identical radiation data, but the lower is overlaid with lines
denoting resonance times for 12 − 14ω (white) and > 14ω (magenta) resonances (there are no
excluded resonances with a smaller or larger number of photons involved for these parameters) for2s − 6p.

109

3.9

Summary
To summarize, we have investigated the near-threshold region of HHG and identified resonance-enhanced radiation associated with states shifting during the pulse. For one particular case – hydrogen driven by an 800 nm, 20 o.c. sine-squared pulse of 3.8 × 1013 W/cm2 – we have thoroughly
demonstrated the process of identifying the responsible resonance (9ω → 3p) and three associateddelays between the time when 3p shifts into resonance and the appearance of radiation. The timingof this radiation remains consistent with predictions within an appreciable variation within theparameter space. In other regimes – including more intense driving lasers, weak mid-IR pulses, andhelium driven by 400 nm pulses – we have identified similar structures, demonstrating the generalimportance of this effect.
While certain features show agreement with an explanation based on semi-classical trajectories, others cannot be thus explained. In particular, we demonstrate a lack of expected CEP
dependence – present at higher intensities and for shorter pulses – and spectral signatures which
do not match the predicted shape. Through multi-parameter scans, in which the resonance time
is kept fixed, we explore the consistency of the 9ω → 3p resonance enhancement. We hypothesizethat breakdown of resonance predictions are associated with low values of UωP (/ 1) and UIPP (/ 0.1),
and that the pulse duration Nτ may play an important role in the timing of radiation (beyond its
impact on the resonance condition).

3.9.1

Future prospects

Resonantly-enhanced near-threshold HHG offers a rich field for further investigation. In this
Chapter, we have focused our attention on a particular feature – enhancement via the 9ω → 3p
resonance – and noted that radiation occurs at multiple distinct delays after the resonance, unlike
for other regimes investigated. If similar multiply-delayed resonance features could be identified(e.g., in helium or with mid-IR pulses), comparison could provide answers to some remaining
questions:110• Do multiple delays appear at consistent delay intervals? (e.g., 0.5 + 1.25n o.c.)• What parameters determine when multiple delays will generate significant radiation, versus
the seemingly more common case of a single delay?
• Do falling-side resonances ever produce significant radiation at second and greater delays?
(For this, we recommend investigation with longer pulses.)
There is of course, another crucial remaining question: If, as we suggest, no trajectory model
is applicable at this and lower intensities, what is the mechanism of radiation? What determines
the delay between resonance and radiation? Are signatures of resonance enhancement present (butperhaps less dominant) in other domains, such as:• radiation at higher harmonics,• from pulses with higher driving intensities, perhaps including cases where tunnel ionization
is expected to play a larger role,• in more complex systems such as molecules or atoms where multiple orbitals contribute(e.g., the 2s and 2p states of neon).
Some tools of investigation that we did not significantly pursue, but that may hold promise:• driving lasers with more (and more continuous) variation of envelope, such as flat-top or
supersine envelopes, such that resonance delays may place radiation within flatter regionsof the envelope,• “spectral profile” plots where radiation is windowed to some particular region for parameterscans,• Floquet analysis or use of numerical basis codes, to investigate deviation from the UP (t)approximation for state shifting,• generalized semi-classical trajectory models (cf. [41]) with multiphoton ionization conditions considering the actual predicted photon energies or states,

111• tuning of state energies separate from ionization potential, through manual modification ofpotential parameters (see Sec. 2.2.5),• and simulations of the response when starting in an excited or superposition state.

Chapter 4Interpolation and Macroscopic Response

We develop a macroscopic description of high-order harmonic (HHG) radiation resulting
from the interaction of atomic systems with an intense laser pulse using ab-initio solutions of
the time-dependent Schrödinger equation (TDSE). After providing some background in Sec. 4.1, inSec. 4.2 we show that for this highly nonlinear process, interpolation can been performed across laser
intensity for a given wavelength, limiting the number of full time-dependent Schrödinger equation
calculations to about one hundred. The impact of the phase of the driving laser on radiation
is explored on the single-atom level in Sec. 4.2.2 and for macroscopic targets in Sec. 4.5. We
summarize the macroscopic harmonic signal computation, which is based on the approach in [106],
in Sec. 4.3. We consider the application of this efficient method for distributions of hydrogen atoms
on a lattice grid (Sec. 4.3.1) and in a gas jet (Sec. 4.4) to demonstrate the successful applicationof the TDSE interpolation method to the thin medium or low gas density regimes that are freefrom longitudinal phase-matching effects. By comparing the results of the present method withtheoretical predictions and against previous results, it is shown that the method based on theinterpolation of ab-initio TDSE solutions leads to reliable macroscopic results. The significantly
reduced computational time as compared to more sophisticated methods opens a path toward theextension of macroscopic high harmonic calculations based on ab-initio microscopic results to more
complex targets and interactions. In Sec. 4.4.2, we demonstrate this via investigations of the nearthreshold regime of the spectra that show that the degree of coherence of the off-harmonic radiation
generated during the pulse is much lower than that of the harmonics. On the other hand, the weak

113
non-harmonic radiation extends to larger divergence angles than the harmonic signals providing
the option for separation of the different signals in this part of the spectrum. Connecting to Ch. 3,
we show in Sec. 4.4.3 that the resonance-enhanced radiation is in fact the portion of off-harmonic
radiation that survives macroscopic propagation, though at relatively weaker intensity than theharmonic radiation. The general applicability of this method is shown for a variety of driving
wavelengths and targets in Sec. 4.6.

4.1

Background
Proper theoretical analysis of HHG experiments requires consideration of macroscopic as-pects, i.e., the coherent build-up of the (microscopic) high harmonic signals from many atoms in
the generating medium. It has been shown (e.g., [14, 107]) that efficient harmonic emission depends
on several macroscopic parameters of the experiment, e.g., beam shape, gas pressure, etc.. The
propagation of the emitted radiation in the medium is described by Maxwell’s equations. Thus,
an exact description of the HHG process requires solution of Maxwell’s equations coupled with thesolution of the TDSE for each of the atoms at the microscopic level. Due to the large numberof atoms in the medium and, hence, the large number of TDSEs to be solved, this task requires
formidable computer resources (see [108–111] for further discussion). In view of the computational
demands, the application range of these methods is very limited and theoretical analysis of most
HHG experiments cannot be realized at this level.
Therefore, a number of models have been proposed and applied in the past to reduce thecomplexity of the problem (e.g., [106, 109–126]) and enable a broader investigation of the parameter
space, incl. by research groups with restricted computational resources. One of these efficient approaches is based on computation of the macroscopic harmonic yield as a superposition of the fieldsof pointlike emitters in the medium via the discrete dipole approximation (DDA), instead of thenumerical solution of the wave equation [106]. Along with an extension of the strong-field approximation (SFA+) to obtain the single-atom high harmonic spectra [53], this SFA-DDA approximationof the macroscopic harmonic build-up has been successfully applied to match experimental data

114
and analyze various aspects of high-harmonic generation, such as generation of bright harmonics inthe keV X-ray regime [31], prediction and application of attosecond vortices from HHG [37, 127],
as well as generation and application of circularly polarized harmonics [38, 103, 128]. The use of
the strong-field approximation or other approximate theories at the single-atom level allows foran efficient calculation of the macroscopic harmonic yield. Although the SFA+ model provides
accurate predictions for the plateau and cut-off region of the spectrum, the application of theseapproximative methods on the microscopic level are limited. For example, they cannot be applied
for an analysis of the signal in the below- and near-threshold region of the high harmonic spectrum.
This part of the spectrum, which has gained an upsurge in interest recently (e.g., [39, 95–101, 129]),
is influenced by the excited state structure of the target atom, which is usually not well reproduced
in approximate single-atom theories. We note that excited states play a role in the control andresonant enhancement of harmonics (see Ch. 3 and [130–135]). The limitation of approximative
microscopic calculations is perhaps even more obvious with more complex targets, like molecules
or nanosystems, which have recently attracted a lot of attention in experiment.
There is a gap in the theoretical analysis of HHG between full ab-initio Maxwell-Schrödinger
simulations, which nowadays cannot be applied to most of the experimental conditions, and thosemethods based on approximative calculations, which are often limited in the accurate description of
microscopic interactions. It is therefore important to develop macroscopic HHG calculations which
include ab-initio microscopic calculations in an efficient way. Such methods will enable extension ofmacroscopic HHG studies to the analysis of the impact of more complex interactions on the spectra
as well as to molecules and nanosystems.
As mentioned above, the impact of excitation on the generated near-threshold part of thespectrum is one of these more complex interactions. A TDSE result for the radiation spectrum
generated in atomic hydrogen is shown in Fig. 4.1. The field-free ionization threshold corresponds
to 8.77ω, where ω is the frequency of the laser pulse, and one can see the typical rich structurein the below- and near-threshold region of the spectrum. In the numerical TDSE solutions at thesingle-atom level, the strength of the signal in between the harmonics – which we denote as off-

115

Figure 4.1: Example of single-atom high harmonic spectrum generated in atomic hydrogen, obtained via numerical solution of the time-dependent Schrödinger equation, exhibiting a rich structure in the near-threshold region of the spectrum between harmonics 5 and 12. The hydrogenionization threshold of 0.5 a.u. = 8.77ω, where ω = 0.057 a.u. corresponds to the central driverwavelength of 800 nm. Other laser parameters: sin2 pulse of 20 optical cycles with peak intensity of
7 × 1013 W/cm2 . Note that the intensity of the spectrum (as well as other spectra in this Chapter)is scaled by the signal at the fundamental frequency.

116
harmonic radiation – appears to be as strong as the harmonic signals. Meanwhile, in experiment,
off-harmonic radiation is usually much suppressed, typically showing up as weak shoulders at theside of a strong harmonic emission (e.g., [96]). Since the radiation is generated during the pulse
and may therefore contain relevant information about ultrashort laser-induced electron dynamics
in excited states, it is interesting to analyze at what strength it occurs in the macroscopic response
and if it can be separated from the harmonic signals.
In this Chapter we present a macroscopic calculation that includes complex effects such as
excitation via TDSE calculations (or other ab-initio methods). To this end, we determine themacroscopic harmonic yield as a superposition of fields of point-like emitters using TDSE results atthe single-atom level. To keep the overall calculation efficient, we consider interpolation of ab-initio
TDSE calculations across intensity. Having in mind that HHG is a highly nonlinear process as afunction of intensity, we thoroughly show that a simple interpolation is indeed sufficient to keep
the number of actual TDSE calculations to a computationally feasible range of about one hundred
for the performance of a macroscopic simulation. Additionally, we demonstrate the accuracy ofa simple approximation of the carrier-envelope phase effect on HHG results. The reliability ofthe macroscopic results is demonstrated by comparison with theoretical predictions for a latticegrid and previous results for the near-threshold regime of HHG in a gas jet under loose focusing
conditions.

4.2

Interpolation
We obtain the microscopic HHG spectra as discussed in Sec. 1.4.3.3, evaluating the dipole

acceleration a(t) using the Ehrenfest theorem
a(t) =

∂− V (r) ,
∂z

(4.1)

where V (r) is the atomic interaction potential. For hydrogen, this is simply − 1r and az (t) = − rz3 .
For other atoms, we consider the derivative of the full potential VSAE . The complex harmonic
response a(ω) is then obtained by taking the Fourier transform of the (un-windowed) dipole accel-

117
eration,Z Ta(ω) =

a(t)e−iωt dt.

(4.2)

0

and the harmonic spectrum P (ω) = |a(ω)|2 .
We note that the interpolation scheme and macroscopic calculation discussed below do notdepend on the particular microscopic solver used – all that is needed is the final a(t) or a(ω).
Thus, the interpolation method can be applied with any preferred TDSE solver or any other abinitio method. This opens the door for exploration of molecular and multielectron systems whose
calculations are computationally prohibitive in combination with other macroscopic methods.

4.2.1

Interpolation across intensity

Interpolation of the TDSE solutions can be accomplished either for a(t) in the time domain
or for a(ω) in the frequency domain. Interpolation in the time domain has the advantage that thedipole acceleration a(t) is purely real. Furthermore, any numerical calculation of the HHG spectra
requires the determination of the dipole acceleration from a numerical simulation. For interpolationin the frequency domain, each numerically determined data a(t) must be Fourier transformed before
the interpolation can be performed. On the other hand, the information content of the harmonic
signal is more efficiently stored in the frequency domain, since one may truncate the spectrum
P (ω) a few harmonics after the HHG cut-off. Depending on pulse duration and numerical time
step, this may decrease memory requirements by more than a factor of 10. Additionally, Fouriertransformation is usually computationally more expensive than interpolation, so if results are only
desired in a single domain (e.g., only spectra are considered), it is most efficient to interpolate in
the desired domain.
The procedure for interpolation used is the same in the time and frequency domains. The
dipole acceleration a(t) (or a(ω)) is numerically calculated for a set of peak intensities over a desired
sampling range. For a given point in time (or frequency), the dipole acceleration (or complex
spectrum) is interpolated across intensities using cubic splines with not-a-knot end conditions. The
impact of the end conditions is minimized by ensuring that relevant intensities are well-contained

118
within the sampling range. The full dipole acceleration (complex spectrum) is then reconstructedfrom these slices. If required, transformation into the frequency (time) domain is completed after
interpolation.
In order to evaluate the success of the interpolation procedure, we have numerically solved
the TDSE for a hydrogen atom at a randomly chosen peak intensity of 9.1431 × 1013 W/cm2
driven by a 20-cycle laser pulse at 800 nm. In Fig. 4.2, we show results of the interpolation for (a)high harmonic spectrum P (ω), (b) dipole acceleration a(t), and (c) error in harmonic phase Φ(ω).
All features of both high harmonic spectrum and dipole acceleration in the exact calculation are
reproduced in the interpolated results. The maximum absolute errors for these specific calculations
are below 5 × 10−6 for the dipole acceleration and below 5 × 10−7 for the harmonic yield. In test
calculations we have found similar results for other driver wavelengths, peak intensities, and pulse
durations.
In the present work, we have chosen to evaluate exact TDSE results for a sample set of peakintensities that are equally spaced on a linear scale, though other spacing (e.g., logarithmic or
Chebyshev) may be preferable in some cases. Results of test calculations suggest that Chebyshev
sampling has benefits especially in other integration schemes than Monte Carlo, potentially requiring 2-3 times fewer sample points. On the other hand, it does limit the re-usability of the sample
calculations for other macroscopic peak intensities.
In Fig. 4.3, we present the maximum absolute error in high harmonic spectra obtained via
interpolation for two different peak intensities (other parameters as in Fig. 4.2) as a function of
the sample spacing dI over the range from 1 × 1012 W/cm2 to 2 × 1013 W/cm2 . As expected for a
cubic interpolation, the error roughly scales as O(dI 4 ) for sufficiently small dI. For larger dI, theerror varies more strongly, and the error is largest when the nearest sampled intensity is far (morethan 2 × 1012 W/cm2 ) from the target intensity. These features are consistent for all wavelengths
tested in the range of 270 − 1600 nm. This result is of particular interest for the long wavelength
regime in view of the unfavorable scaling of computation times for HHG calculations with increaseof the wavelength.

119

Figure 4.2: (a) Comparisons of results of interpolation (dashed lines) and exact numerical results(solid lines) for high harmonic spectrum and (b) for dipole acceleration, and (c) error in harmonic
phase, generated by a sin2 pulse of 20 optical cycles at wavelength of 800 nm and peak intensity of
9.1431×1013 W/cm2 . Interpolation performed on data with intensities of 0.1, 0.2, 0.3, ..., 10.0×1013
W/cm2 . ((a) taken from [6].)120

Figure 4.3: Trend of maximum absolute error in interpolated high harmonic spectrum (as compared to numerical TDSE result) as a function of the intensity sampling spacing dI (i.e., sampledintensities dI, 2dI, . . . , ∼ 1014 W/cm2 ). Shown are exemplary results at peak intensities 3.97 × 1013
W/cm2 (blue circles with dashed line) and 9.14 × 1013 W/cm2 (red circles with dotted line). Other
laser parameters as in Fig. 4.2. Note that at sufficiently small intensity spacings, the error roughly
drops as O(dI 4 ) (black dashed lines).

121
While the choice of intensity sampling spacing depends on the desired accuracy, for 800 nm
we have found dI = 1012 W/cm2 to be fully sufficient. In general, from our test calculations atother wavelengths, we have found that the spacing required relates roughly to the maximum extent
of the harmonic spectra, i.e., 3.17UP /ω ∝ I/ω 3 . Thus, as a rule of thumb the sampling spacing forthe cubic interpolation should be (in W/cm2 ): ω 3
800 3
12
dI ≈
× 10 =× 1012
0.057
λ

(4.3)

where ω is in a.u., while λ is in nm.
While our primary focus in this Chapter is to use the interpolation method to facilitate
macroscopic HHG calculations, we note that the success of the interpolation method is interesting
in itself since HHG is a highly nonlinear process and depends strongly on intensity. As we haveshown, cubic interpolation produces quantitatively good agreement with exact spectra at moderateintensity sampling (and qualitatively indistinguishable spectra can be generated with coarser sampling). While perhaps not surprising from the numeric side, the implication that the results of sucha highly nonlinear process can be distilled to order-100 samples for description of the whole range of
relevant intensities (with saturation of ionization well before 1015 W/cm2 ) is not a-priori obvious.
In view of the success, it is likely that a similar technique can be applied for the dependence of
HHG spectra on other parameters as well.

4.2.2

Variation of CEP

In addition to different driving intensities, it is important to consider the effects of variationof the carrier-to-envelope phase (CEP), φ, on the spectrum. Especially in short pulses, CEP is aparameter that controls many strong-field applications such as isolated attosecond pulse generation[136], laser-driven electron dynamics in atoms, molecules, and nanostructures [137–140], or generation of relativistic electron beams [141]. For long pulses, the impact on the microscopic spectral
intensity is minimal (Fig. 4.4), but the spectral phase effects cannot be discounted in macroscopic
summation.

122
Before proceeding, we provide a note on interpreting phase plots. Plots of phase tend toinclude large spikes that have little physical significance. The first reason is that phase has little
meaning when the magnitude is small: +10−4 and −10−4 are much the same when features ofinterest are on the scale of ±1. Secondly, phase is defined modulo 2π, and smoothing out 2π
jumps is not always straightforward or useful. The simplest choice of “unwrapping” phase byminimizing the difference between adjacent points leads to spectral phases spanning hundreds ofradians, obscuring most features. On the other hand, restricting phases to [0, 2π) can appear to
break larger trends. For this reason, we select CEPs within a specific range that allows simpler
presentation of spectral phase without manual addition of 2π’s. The results are consistent for larger
CEP.
In SFA, the spectral phase is given by:Φ(ω, φ) = Φ(ω, φ = 0) + φH(ω)(4.4)

where H(ω) is the harmonic number rounded to the nearest odd integer [106]. In Fig. 4.5, we show
exact TDSE results for the phase difference∆Φ/φ = Φ(ω, φ) − Φ(ω, φ = 0) /φ(4.5)

compared to the SFA prediction ∆Φ/φ = H(ω). For both intensities, we see a factor between CEPand spectral phase which does not depend on CEP (i.e., the scaled lines in Fig. 4.5 overlap):∆Φ/φ = f (ω)(4.6)

for some function f (x). For the lower intensity, there is good agreement with fSF A (ω) = H(ω).
For the higher intensity, H(ω) approximates the true dependence, but the specifics are more complicated.
While not further investigated in this thesis, we do present results for a short (4 o.c.) pulse
at the lower intensity (Fig. 4.6). The spectral intensity differences are larger, but mostly not
impacting the harmonics themselves. The spectral phase varies more dramatically (especially inthe near-threshold 6 − 10ω region of interest), and in places does depend on the particular CEP in123

Figure 4.4: Spectra for a selection of sample CEPs, generated by a sin2 pulse of 20 optical cycles
at wavelength of 800 nm and peak intensity of (top) 3.8 × 1013 W/cm2 and (bottom) 7.0 × 1013
W/cm2 .

124

Figure 4.5: Spectral phase difference for a selection of sample CEPs, generated by a sin2 pulse of 20
optical cycles at wavelength of 800 nm and peak intensity of (top) 3.8 × 1013 W/cm2 and (bottom)
7.0 × 1013 W/cm2 . (Top panel taken from [6].)125
a way not seen for longer pulses. Still, for the majority of the spectrum, the CEP effects are similar
to those for longer pulses – the radiation intensity is unaffected and the SFA approximation to thephase is good. We note parenthetically that the driving laser field was not frequency-corrected [85],leading to the harmonic blue-shift.
In Sec. 4.5, we will investigate the impact of the spatial CEP variation on the macroscopic
radiation using the ∆Φ/φ = H(ω) approximation. With the addition of Gouy phase, the effects
are small and smooth enough that we predict the true dependence would not significantly affectresults for long pulses. For highly-focused pulses where an additional focal phase [142, 143] becomessignificant (see Sec. 4.5), we expect the same approximation to hold, but recommend future study
involve consideration of interpolation across CEP. More careful treatment will also be needed forcircularly and elliptically polarized laser pulses, as well as when multiple pulses contribute to HHG.

4.3

Macroscopic summation
To determine the macroscopic radiation signal at a far-field detector, we consider the thinmedium or low gas density regimes that are free from longitudinal phase-matching effects. We
define the arrangement and coordinates as follows: The driving laser propagates in the ẑ direction,
has linear polarization ŷ, and we consider only radiation polarized in this direction. The far-field
detector is an arc in the xz-plane with its position defined by the angle θ from ẑ (propagation)
to x̂. We then follow the approach used in [106], in which the macroscopic yield is obtained as
the superposition of the fields generated at different points in the medium. This approximation
to the full Maxwell solution relies on the dipole approximation and the assumption that generated
radiation does not interact with the medium.
The spectral distribution of the total radiation generated by a number of atoms, located atrj (j = 1, 2, 3, ...), is given by [106]:E(rd , ω) =ω
1 −i ω |rd | Xe caj (ω)e−i c [rj ·(r̂d −ẑ)] ŷ ,2
c |rd |

(4.7)

jwhere aj (ω) is the dipole acceleration in frequency domain of a single atom responding to the driving

126

Figure 4.6: (top) Spectrum and (bottom) spectral phase difference for a selection of sample CEPs,
generated by a short sin2 pulse of 4 optical cycles at wavelength of 800 nm and peak intensity of
3.8 × 1013 W/cm2 .

127
laser pulse. Since this is calculated without considering the delay of the driving laser reaching theradiator, (rj · ẑ)/c, the corresponding additional phase is added. Furthermore, it has been assumed
that the relative location of the detector, denoted by rd , is far away from the individual atoms, i.e.,1|rj |  |rd | and, hence, |rd − rj | ≈ |rd | − rj · r̂d and |rd −r≈ |r1d | . Since the prefactor in Eq. (4.7)
j|universally scales the results, we drop it from the computations. As usual, radiation is reported as
spectral intensities, |E(rd , ω)|2 .
We approximately include the effect of the Gouy phase ζ(z) = − atan(z/z0 ) [144], on Hthharmonic radiation as an additional phase H(ω)ζ(z), where z0 is the Rayleigh range, as in [106]. As
shown in the previous section, for 20 o.c. pulses this approximation holds well. Additionally, for theloose focusing geometry used below, ζ ranges up to ≈ ±0.2 rad within the gas jet, and the effect onnear-threshold harmonics – the main focus of the applications below – is minimal. Specifically, as
seen in Sec. 4.2.2, the CEP effect is mostly restricted to a consistent additional phase scaling withharmonic number. This has the effect of adding a phase mismatch for higher harmonics, leading
to suppression. Centering the gas jet off the laser focus allows for better phase matching when
the Gouy and intrinsic phases compensate [106]. In Sec. 4.5, we discuss more closely the effects ofGouy phase as well as another phase term related to highly-focused laser pulses.
Using rd = rd (cos(θd )ẑ + sin(θd )x̂) and rj = rj (cos(θj )ẑ + sin(θj )x̂) +yj ŷ, Eq. (4.7) can thenalso be written as:
E(θd , ω) =Xωaj (ω)e−i c rj [cos(θj )−cos(θd −θj )] ŷ.

(4.8)

j

The sum in Eq. (4.7) or (4.8) can be efficiently performed if the number of atoms does not become
too large. However, in a typical high-harmonic gas jet experiment, even at low-pressure densities,
the number of atoms is typically of the order of 1012 or larger. To simplify the computation in thesecases, it has been proposed [106] to perform the sum using Monte Carlo sampling, and further tosplit the interaction volume into discrete cells, each containing a number of radiating atoms. The
atoms in the cell approximately interact with the same external field if the number of atoms in thecell is large enough to approximate their density by a continuous distribution and the size of the128cells is small enough to approximately have a uniform intensity distribution within the cell. The
total field at the detector is then given by [106]:E(θd , ω) =

X

ω

nj

aj (ω)e−i c rj ·[(cos(θd )−1)ẑ+sin(θd )x̂]×e− 2 c2 σ (1−cos(θd )) ŷ ,

j
1 ω2

2

(4.9)

where the sum goes over the number of cells in the interaction volume, nj = ρ(rj ) is the gas
density at the center of the cell and the last term is a form factor that takes into account intra-cellinterference assuming a Gaussian charge distribution in the cell. Note that the use of cells does not
impact the radiation detected along the laser axis (θd = 0).

4.3.1

Test case: Lattice grid

As a first test, we consider a set of hydrogen atoms placed in a regular lattice grid and verify
that the angle-resolved far-field radiation spectra exhibit the expected interference pattern. To thisend, we apply the TDSE interpolation method to a simple test model in which hydrogen atoms,
initially in the ground state, are distributed uniformly on a lattice with a separation of 300 nm(typical atomic spacings in optical lattices are in the range of 256 nm - 700 nm [145, 146]). We haveconsidered a driving laser with a central wavelength of 800 nm (ω = 0.057 a.u.), a peak intensity of
9×1013 W/cm2 , and a sin2 envelope of 20 cycles full width. The spatial profile of the laser is chosen
to be a Gaussian beam with a beam waist of 30 µm (Rayleigh range z0 ≈ 3500 µm). One hundred
separate TDSE calculations were performed at peak intensities in the range 0.01 × 1014 W/cm2 to
1.00 × 1014 W/cm2 in steps of 1012 W/cm2 and the interpolation of the results is completed in thefrequency domain to reduce memory requirements. Given the small number of radiators here, we
sum over all atoms, without using Monte Carlo sampling or finite cells.
Considering a regular lattice geometry with atoms arranged along and perpendicular to the129

Figure 4.7: Angle-resolved macroscopic high harmonic spectra from uniformly-separated hydrogenatoms in different configurations: (a) row of 100 atoms along the propagation direction, (b) row of
100 atoms perpendicular to the propagation direction, and (c) lattice grid of 100 × 100 atoms in
the plane spanned by the combination of (a) and (b). A Gaussian beam (beam waist 30 µm) with
20-cycle duration at 800 nm and 9 × 1013 W/cm2 and an atom separation of L = 300 nm wereconsidered. The red and white lines represent the constructive interference conditions for atoms
located along and perpendicular to the propagation direction, respectively (see Eqs. (4.11) and(4.10)).

130

Figure 4.8: Comparison of spectra taken from the results in Fig. 4.7(c) along θd = 0.65 rad (bluesolid line) and in forward direction (θd = 0 rad, black dashed line).

131
laser propagation direction, radiation from the atoms is expected to add up constructively if
θd
θd



λem
= arcsin n
and
L
!
r
m λem
= 2 arcsin,
2 L(4.10)(4.11)

where n and m are integer numbers, L denotes the separation of the atoms and λem is the wavelengthof the emitted radiation. The first (second) condition arises from the interference in emission fromatoms perpendicular to (along) the laser polarization direction, and can be derived from Eq. (4.8).
Thus, we expect the process of high-order harmonic generation from uniformly spaced,independently-responding atoms to produce a particular pattern in the angle-resolved spectrum.
To test the TDSE interpolation method, we have performed macroscopic calculations for 104 hydrogen atoms on a 100 × 100 lattice. The results in Fig. 4.7(c) show the complete spectrum inthe forward direction (at θd = 0) and a pattern of bright spots at other angles. As expected fromthe constructive interference conditions, the number of bright spots in a given detector angle range
increases as the (harmonic) frequency of the emitted radiation increases (i.e., the wavelength λem
decreases) and the spectra show the expected symmetries about θd = 0, ±π/2. Indeed, the brightspots occur at the vertices of curves representing the two interference conditions. To illustrate
this, we show in the other panels of Fig. 4.7 the spectra resulting from configurations in which 100
atoms are positioned in a row either along (panel (a)) or perpendicular to (panel (b)) the laser
propagation direction. Furthermore, the curves representing the respective interference conditions
are plotted in each spectrum. It is seen that the numerical results indeed agree with the predictions
based on the constructive interference conditions.
While these calculations primarily serve as a test case for the macroscopic TDSE interpolation
method, we may note that optical lattices are used as platforms in ultracold atom physics to simulatequantum many-body problems (for reviews, see [147–151]). Recently, ultracold atoms have been
exposed to the field of an ultrashort laser pulse to observe the response of 87 Rb atoms to strong
fields [152], which may initiate more combined research in these fields in the future. Here, we notethat the spatial separation of the bright spots in the angle-resolved spectrum generated from atoms

132
in a lattice potentially allows for the selection of radiation at specific higher harmonics. This is
illustrated in Fig. 4.8 in which we compare the spectrum at θd = 0.65 rad (blue solid line) withthat in the forward direction (black dashed line). In the off-axis spectrum, only the fundamental
exceeds 0.1% of the 13th harmonic’s intensity while the rest of the spectrum is strongly suppressed.

4.4

Atomic hydrogen gas jet
In a common experimental HHG set-up, a large number of atoms in a gas jet contribute to

the macroscopic signal. For this section, we therefore consider a gas jet of hydrogen atoms witha density of 1018 cm−3 . The gas jet is modeled as a cylinder aligned perpendicular to the laser
propagation; the width is given by a Gaussian distribution with σ = 500 µm and is centered at thelaser focus. We use the strategy to split the interaction volume in cells and utilize Monte Carlo
methods for the sampling of the cell locations, discarding locations where the intensity is below agiven cut-off (here, I(rj ) < 0.032Ipeak ) as the corresponding generated radiation does not contributesignificantly to the total radiation. We have considered macroscopic high-harmonic generation fromhydrogen atoms generated by a 20-cycle sin2 pulse at 800 nm, with a variety of peak intensities
up to 9 × 1013 W/cm2 . As above, the spatial profile of the laser is chosen to be a Gaussian beamwith a beam waist of 30 µm (z0 ≈ 3500 µm). One hundred microscopic TDSE calculations were
performed at intensities in the range 0.01 × 1014 W/cm2 to 1 × 1014 W/cm2 in steps of 1012 W/cm2 .
We note the need for single-atom calculations above the highest intensity of the macroscopic lasers
considered, to allow for accurate interpolation.

4.4.1

Convergence

For the case of a gas jet, we find similar convergence behavior to the SFA+ calculations inRef. [106]: the off-axis harmonic radiation appears to converge faster for larger cells (Fig. 4.9),
and the results appear to be converged for 105 radiators regardless of cell size (Fig. 4.10). In
general, we observe that good on-axis convergence requires a larger number of cells than required
for angular convergence, at least for the near-threshold harmonics of interest in this thesis. Thus,133
cell size – which has no effect for on-axis radiation – does not noticeably impact the speed of overall
convergence, and we use point-like radiators for all calculations. We note that the convergence in thecalculations of the response using 105 cells has been achieved from 100 actual TDSE calculations at
sample intensities. Since the macroscopic calculations are relatively computationally inexpensive,
results of this Chapter consist of calculations with 5 × 105 atomic radiators unless specified.
We must emphasize that the required number of cells is not universal, and may vary greatly
for other lasers or targets. We have found that both interaction of UV laser light with helium
and mid-IR lasers with hydrogen require on the order of 107 cells to achieve the same levels of
convergence as above (see Sec. 4.6). More rapid or complex CEP variation also requires morecells for convergence (Sec. 4.5). Preliminary studies suggest alternative cell placement (e.g., quasirandom sampling) may improve convergence somewhat.

4.4.2

Near-threshold harmonics

As discussed at the outset of this Chapter, approximate single-atom methods, such as thestrong-field approximation and its extensions, provide good quantitative predictions for the highorder harmonics. On the other hand, in the region near and below the ionization threshold, singleatom TDSE results and experimental data often include a structure besides the harmonic peaks,
which often cannot be reproduced by approximate methods due to the incomplete description of
the atomic excited state structure.
In Fig. 4.11 we compare results of (a) microscopic single-atom TDSE and (b) macroscopic
TDSE interpolation calculations as a function of peak intensity for on-axis emission. In the results it
is seen that at the single-atom level there is strong radiation emission at near-threshold frequencies
in between the harmonics (the field-free ionization energy is at 8.77ω). The relative strength of theharmonic and off-harmonic radiation emission strongly fluctuates as a function of peak intensity. As
examples, we present spectra at two different peak intensities ((c) 3.8 × 1013 W/cm2 , (d) 7.0 × 1013
W/cm2 ) on the right of Fig. 4.11. While at the higher intensity (panel (d)) there is a strong singleatom 9th harmonic emission (black dashed line), the emission at the same harmonic is suppressed

134

Figure 4.9: Angle-resolved macroscopic radiation spectrum for a gas jet of hydrogen atoms driven
by an 800 nm laser with peak intensity 7.0 × 1013 W/cm2 for different numbers of cells (rows) andcell sizes (columns): (top) 103 cells, (middle) 104 cells, and (bottom) 105 cells. (left) point-likeradiators, (middle) 3 µm cells, and (right) 5 µm cells.

135

Figure 4.10: Convergence of macroscopic on-axis radiation with number of cells (cell width has noeffect), generated by a sin2 pulse of 20 optical cycles at wavelength of 800 nm and peak intensity of(top) 3.8×1013 W/cm2 and (bottom) 7.0×1013 W/cm2 . Vertical dashed lines denote the ionization
potential (left line) and HHG cut-off (right line).

136

Figure 4.11: Comparison of single-atom TDSE results and macroscopic TDSE interpolation resultsfor high harmonic generation in hydrogen atoms. On the left, the (a) single-atom and (b) macroscopic on-axis spectrum is shown as a function of peak intensity. The spectra on the right show
the comparison (single atom: black dashed line, macroscopic: red solid line) at peak intensities (c)I0 = 3.8 × 1013 W/cm2 and (d) I0 = 7.0 × 1013 W/cm2 .

137(as compared to radiation emitted at frequencies just below the 9th harmonic) at the lower intensity(panel (c)). In general, in the near-threshold regime, the magnitude of the off-harmonic radiation
yield is comparable to that of the harmonic radiation.
In contrast, the macroscopic results (panel b) show emission spectra with a much clearer
structure; the harmonics are strongly enhanced as compared to radiation emission at other frequencies at all intensities. This is also evident in the macroscopic responses at the two selectedpeak intensities. The radiation signal in between the harmonics is suppressed by 1-2 orders of magnitude as compared to the harmonic emission. For example, the feature close to the 9th harmonic
appears as a shoulder to the harmonic. Thus, in the macroscopic spectra the high harmonic peaksin the near-threshold region become as prominent as those in the rest of the spectrum, throughout
the plateau to the cut-off, similar to typical observations of high-harmonic spectra in experiment(see, e.g., [96]).
The origin of the suppression of the off-harmonic radiation is related to the differences in
the coherent build-up of the distinct parts of the spectrum. This can be seen from Fig. 4.12 inwhich we present a comparison of calculations where the interpolated single-atom TDSE results
are either added up coherently (solid red line, as Eq. (4.8)) or incoherently (dashed blue line, withP
P (ω) = j |aj (ω)|2 ). It can be clearly seen that it is the lack of coherence of the radiation at
non-harmonic frequencies that is the primary reason for the strong suppression in the macroscopic
result.
Although in the on-axis spectra the off-harmonic near-threshold radiation is strongly suppressed, this part of the emitted radiation can be separated from the harmonics since it has adifferent angular character than the harmonic radiation (cf. Fig. 4.13). The high-harmonic lines
show the expected on-axis emission, which becomes more focused in the forward direction for higher
harmonic numbers. In contrast, the off-harmonic emission extends to larger divergence angles and,around some harmonics, forms halo-like patterns in the angle-resolved macroscopic spectrum. Suchspatial-spectral characteristics of off-harmonic structures in the low- and high-energy part of thespectrum are in agreement with previous observations [95, 96, 153–157], showing that the present

138

Figure 4.12: Comparison of macroscopic calculations in which the single-atom TDSE results areadded coherently (solid red line) and incoherently (dot-dashed blue line) for high-order harmonic
generation at I0 = 3.8 × 1013 W/cm2 .

Figure 4.13: Angle-resolved near-threshold part of the macroscopic radiation spectrum for a gas
jet of hydrogen atoms driven by an 800 nm laser with peak intensity 3.8 × 1013 W/cm2 .

139
efficient interpolation method leads to similar results and conclusions as previous sometimes moresophisticated simulations.
In Fig. 4.14, we show wavelet analysis (see Sec. 1.5.3) for the radiation (a) on-axis and (b)at 5 mrad off-axis. The on-axis radiation is dominated by consecutive odd harmonics, which – asexpected – are emitted at times symmetrically about the center of the pulse. On the other hand, theoff-axis wavelet, which primarily contains the non-harmonic emission, shows a different and more
complex time dependence. In particular, we note that several off-harmonic features appear firstnear the peak of the pulse and shift downward in energy at the trailing edge. The curvature of theshift roughly scales with the field-shifted energy of the excited states, and the wavelet results thus
suggest that the off-axis radiation in the near-threshold part of the radiation spectra are related to
the resonance features, discussed in Ch. 3.

4.4.3

Appearance of resonance features in macroscopic response

We note that the results presented in Figs. 4.11 to 4.13 have been obtained for propagationof the wavefunction until the end of the driving laser pulse. Thus, the radiation discussed above,
including the off-harmonic spectral features, has been generated during the interaction between
the laser pulse and the hydrogen atom. It is therefore unlikely that the structures can be related
to free induction decay following coherent excitation in the atom. Instead, as discussed in Ch. 3,
we understand these features to be the result of field-shifted resonances enhancing HHG emission.
By inverse Fourier transforming the macroscopic E(ω), we are able to perform the same sorts ofanalysis for macroscopic results as seen in Ch. 3 for microscopic results. Given the interpolation
requirement of scanning intensity for any given laser parameters, it is not currently feasible toperform macroscopic scans across other parameters (e.g., λ, Nτ , Z). In the future, interpolation
across multiple laser parameters may enable efficient exploration of the broader space.
In Fig. 4.15, results for radiation near 9ω as a function of intensity are given for microscopic(as Fig. 3.10) and macroscopic calculations. While the majority of off-harmonic radiation is strongly
suppressed, some features survive macroscopic summation. In particular, macroscopic spectra near

140

Figure 4.14: Wavelet analysis of the macroscopic radiation for a gas jet of hydrogen atoms driven by
an 800 nm laser with peak intensity 7 × 1013 W/cm2 at detection angles (a) θ = 0 mrad (on-axis),and (b) θ = 5 mrad. The shifted energy of the 2p, 3p, and 4p states relative to the 1s (ground)state are denoted by red lines, while the region containing the higher Rydberg states (shifted 5p tothreshold) are shaded in red. Wavelet parameter Ω0 = 54.

141
2.0×1013 W/cm2 display similar radiation peaks splitting off 9ω and broad radiation in the region offield-free Rydberg energies, as seen in microscopic results. While the intensity of these macroscopic
features is less than the on-harmonic radiation, it is still significant.
In Fig. 4.16, we consider temporal profiles as outlined in Sec. 1.5.1, again comparing microscopic and macroscopic radiation. The macroscopic signatures of resonantly-enhanced harmonic
emission are quite prominent, with a few key differences. Firstly, radiation just after the peak of
the pulse is consistently of similar strength to the strongest distinguishable resonance features. Second, resonantly-enhanced radiation occurs for intensities somewhat above the intensities predicted.
Both features are explained by the range of intensities contributing to any given macroscopic result.
“Normal HHG” occurring just after the peak appears consistently for wider intensity ranges thandoes resonantly-enhanced emission (since resonances move). And, only a small minority of atoms
interact with the peak intensity; the signal is thus dominated by atoms interacting with a slightly
weaker field. This also leads to a broadening (in time and driving intensity) of resonance features– for a given peak intensity, different portions of the gas jet will shift into resonance at different
relative times. As a result of resonant enhancement, radiation near 9ω spans varying temporal
ranges according to the driving intensity. Therefore, understanding resonant effects will be crucial
to sculpting UV pulses.

4.5

Effects of macroscopic phase1
Nonlinear processes in ultrafast light-matter interaction are induced at high laser peak in-tensities and thus require the application of focused short-pulsed lasers with a broad frequency
bandwidth. In contrast to above assumptions of a monochromatic Gaussian beam which experiences a Gouy π-phase shift across the focus [144], broadband pulses have a different spatial phasedependence. This more general phase distribution has been derived and studied by Porras andco-workers [142, 143]. The exact distribution depends on the so-called Porras factor, which needs
to be determined for a given laser system. It creates a spatially-dependent CEP. The impact of1

Some results of this Section are also presented in a submitted manuscript [6].

142

Figure 4.15: Comparison of (top) single-atom results and (bottom) on-axis macroscopic TDSE
interpolation results for near-threshold HHG spectra as a function of peak intensity for an 800 nm,
20 o.c. pulse interacting with atomic hydrogen. Results for each intensity have been normalized
to the maximum signal in the region between 8.75ω and 9.5ω. Horizontal red dashed lines denotefield-free excited state transition energies and vertical dashed lines identify intensities selected inSec. 3.4: 2.0 × 1013 , 3.0 × 1013 , 3.8 × 1013 , 5.6 × 1013 , and 7.0 × 1013 W/cm2 .

143

Figure 4.16: Comparison of (top) single-atom results and (bottom) on-axis macroscopic TDSE
interpolation results for temporal profile of radiation within 1ω of the 9th harmonic as a functionof driving pulse intensity in hydrogen for an 800 nm, 20 o.c. pulse. For each intensity, the results
have been normalized to a maximum of 1 for sake of comparison. A vertical dashed line denotes
the peak of the laser and horizontal dashed lines identify intensities selected in Sec. 3.4: 2.0 × 1013 ,3.0 × 1013 , 3.8 × 1013 , 5.6 × 1013 , and 7.0 × 1013 W/cm2 . Lines are drawn for delayed resonances
as outlined in Table 3.1.

144
the focal phase on electron backscattering at nanoscale metal tips [158] and photoelectron spectra[159] in few-cycle laser pulses has been demonstrated recently.
Since in HHG the coherent addition of radiation controls the results of macroscopic summation (cf. Fig. 4.12), one can expect to observe traces of the phase distribution in the radiation. We
present macroscopic high harmonic generation results that indicate significant differences in the
angular resolution of the emitted harmonics, depending on the spatial phase distribution across thefocus. Specifically, our results show a distinct interference pattern in the angular distribution for
focal phase distributions with negative Porras factors.
The general form of the spatial phase distribution for broadband laser pulses is given by[142, 143]:∆ΦF (z, r) = − arctan

z
z0


2 

r
1 − 2 w(z)


+ g0

z0
z
z0 + z(4.12)

where
g0 =

dZR (ω)ω0dω ω0 ZR (ω0 )(4.13)

is the Porras factor andZR (ω) =ωW 2 (ω)
2c(4.14)

is the frequency-dependent Rayleigh range (z0 ≡ ZR (ω0 )), z and r are the positions along and transq
verse to the propagation of the laser, ω0 is the central frequency of the laser, w(z) = W0 1 + ( zz0 )2 isthe beam radius of the center frequency, W (ω) is the frequency-dependent input waist, W0 = W (ω0 )is the central frequency beam waist at the focus, and c is the speed of light. The first term in
Eq. (4.12) is the Gouy phase for monochromatic beams, which gives a longitudinal phase dependence along the laser propagation direction, while in the transverse direction the Gouy phase isconstant. The second term scales with the Porras factor g0 and describes the difference to the Gouy
phase for beams with a frequency bandwidth.
Fig. 4.17 shows a comparison of the spatial phase profile for (a) g0 = 0 (Gouy phase) and(b) g0 = −2. We have focused on situations with negative Porras factor, which correspond to
the geometry of laser systems based on hollow-core fiber-compressors [159]. The specific value of145

Figure 4.17: Spatial distribution of carrier-envelope phase in a focused broadband Gaussian laser
pulse: (a) Gouy phase (g0 = 0) and (b) focal phase with g0 = −2. For the latter, reference lines
are drawn showing (magenta) position where laser intensity falls below 1% of peak, (cyan) position
where gas density ρ falls below 10% of peak for a gas width of 0.23z0 , and (white) estimates of theeffective slit locations based on interference patterns (0.925 times the transverse node separationnear z = 0 in Eq. (4.12)). (Taken from [6].)146
g0 = −2, chosen in the majority of our studies, has been reported in a recent measurement [158].
For negative values of g0 , the 2nd focal phase term has the same sign as the Gouy phase on-axis,
but changes sign at certain radial distances, which depend on g0 and z. This implies that for the
total focal phase at each off-center position along the propagation distance there are two points in
the radial direction, symmetric over r = 0, with ∆ΦF = 0. Furthermore, these points are in-phase√with all the central points at z = 0 and for the location of the points limz→0 rzero (g0 , z) = W0 / 2.
One may therefore expect that such a focal phase distribution imprints an interference structure
onto the angular distribution of the harmonic signals transverse to the propagation direction.
In Fig. 4.18, we present simulated angle-resolved macroscopic high harmonic spectra generated by a driving laser with a central wavelength of 800 nm (ω0 = 0.057 a.u.), a peak intensityof 7 × 1013 W/cm2 , and a sin2 -envelope of 20 cycles full width. The spatial profile of the laser
has been chosen to be a Gaussian beam with a beam waist of 30 µm (z0 ≈ 3500 µm). We haveassumed a gas jet with a Gaussian density distribution along z, centered about z = 0 and with awidth of σz = 800 µm ≈ 0.23z0 , and constant density along x (the radial direction). Results havebeen obtained assuming focal phase distributions with (a) g0 = 0 (Gouy phase distribution) and(b) g0 = −2. In the present simulations, 107 single-atom results have been used to obtain eachconverged macroscopic HHG spectrum.
The comparison reveals a significant difference in the angular distribution of the belowthreshold harmonic lines for the two phase distributions. The spectrum obtained with negative
Porras factor clearly exhibits the expected interference pattern in the harmonic lines, while there isno such signature in the spectrum calculated for the Gouy phase distribution. The position of theside maxima in the interference pattern can be roughly approximated by the double-slit formula(i.e., Eq. (4.10)):λem
θd (λem ) = arcsin nL,(4.15)

where n is an integer, L denotes the separation of the slits perpendicular to the laser propagation
direction and λem is the wavelength of the emitted radiation. The white lines in Fig. 4.18 show

147

Figure 4.18: Comparison of simulated angle-resolved macroscopic high harmonic spectra using (a)Gouy phase (g0 = 0) and (b)√focal phase with g0 = −2. Lines denote peaks of a double-slit patternwith slit separation of 0.925 2W0 ≈ 39µm in the direction perpendicular to the laser propagation.
Laser and gas jet parameters are given in the text. (Taken from [6].)148
predictions for the side maxima, fitted to match the numerical data using L = 0.925L0 . Thiseffective slit separation is slightly smaller than the distance between the two radial positions L0 =√2W0 for ∆ΦF = 0 near z = 0. We further note that the off-harmonic radiation at energies in
between the harmonics does not show an interference pattern. This is due to the incoherent natureof this radiation, which results from fluorescence due to resonant excitation of the hydrogen atom[160].
A key factor for the occurrence of the interference effect is that for nonzero g0 -values thephase shift equals zero at radial positions in the focus with z 6= 0. This raises the question of how
the interference pattern depends on parameters of the set-up. Results of those studies are presented
in Fig. 4.19. In the plots, we show how the positions of the first minimum (circles) and the firstmaximum (squares) in the interference pattern of the 7th harmonic depend on (a) the width of thecentral gas jet σz , (b) the Porras factor g0 , and (c) the beam waist W0 . Other parameters are kept
the same as in Fig. 4.18.
The positions of the first two extrema remain almost unchanged for variations of the central
gas jet width and the Porras factor. We note that an increase of the gas jet width leads to adecrease in the contrast between the maxima and minima in the interference pattern until foreven larger jet widths, i.e., σz > 0.5z0 , the pattern visually disappears. This is likely due to theenhanced impact of out-of-phase contributions from locations in between the nodal points which
tend to destructively interfere with the contributions from points with ∆ΦF = 0. Within the limits
for the jet widths established in Fig. 4.19(a), the position of the nodal points does not change
significantly along the propagation axis. This explains the stability of the extrema over the range
of jet widths. The Porras factor g0 controls the convexity (concavity) of the curve along which thenodal points appear as a function of z. However, near the center of the focus, the change is not
significant enough to impact the position of the extrema in the interference patterns for σz ≈ 0.23z0(cf. Fig. 4.19(b)). On the other hand, for a variation of the laser beam width W0 (Fig. 4.19(c)), the−1/2

locations of the interference maxima approximately scale with W0, as expected for the simplified

interpretation based on the double-slit interference pattern. Small variations from the trend are149

Figure 4.19: Position of first minimum (circles) and first maximum (squares) in the interference ofthe 7th harmonic as a function of the (a) width of the central gas jet σz , (b) g0 , and (c) W0 . Other
parameters are kept the same as in Fig. 4.18. The line in the bottom panel shows the prediction for
the first side maximum based on the double-slit formula with estimated slit separation of 0.925L0 .
(Taken from [6].)150
seen at the smallest values of W0 . Since z0 ∝ W02 , this can be understood again as likely due toout-of-phase contributions from the relatively wide gas jet.
In Fig. 4.20, we demonstrate the suppression of higher harmonic radiation on-axis as a resultof the macroscopic phase. As before, the single-atom response (black) has relatively equal intensity for harmonics through the cut-off, as well as for the off-harmonic radiation. Including only
intensity variation and geometric phase, the macroscopic result for a (non-physical) laser with nospatial phase variation (purple) strongly suppresses most off-harmonic radiation, and decreases theintensity of the highest harmonics (> 15ω) relative to the plateau. Including the Gouy phase of a
monochromatic Gaussian laser (red), radiation is suppressed proportional to the harmonic order,
but the structure of the off-harmonic radiation does not greatly differ. When the additional focal
phase is included (green), on-axis radiation is further suppressed for higher harmonics, resulting
in radiation at the cut-off 10 orders of magnitude weaker relative to the fundamental than that
present in the single-atom response. When present, this additional focal phase must be consideredin phase-matching.
Within the numerics of our macroscopic simulations, the focal phase is also significant.
Fig. 4.21 demonstrates the increased convergence requirements when including the focal phase(bottom Panel) as compared to calculations with only the Gouy phase (top Panel). The focal
phase requires orders of magnitude more cells for convergence, especially for harmonics above thethreshold. We attribute this to a combination of lower radiation intensity – requiring more cells tobring the noise floor below regions of interest – and more rapid phase variation within the gas jet– requiring more cells to well-sample the phase structure.

4.6

Other investigations
In addition to our core focus of hydrogen driven by an 800 nm laser pulse, we investigatedother wavelengths and targets, with parameters summarized in Table 4.1 (for example chosen
intensities).

151

Figure 4.20: Comparison of single-atom results (black dotted) and macroscopic TDSE interpolation
results for different driving laser CEP profiles: constant phase (purple), Gouy phase (red), focal
phase with g0 = −2 (green). For a gas jet of hydrogen atoms driven by an 800 nm laser with peakintensity 7.0 × 1013 W/cm2 .

Parameter
Target atom
Wavelength (nm)Intensity (W/cm2 )Duration (o.c.)
IP (a.u.)
IP /ω
UP /ω
Keldysh γ

Hydrogen800
3.8 × 1013
20
0.50
8.8
1.5
1.7

Hydrogen1600
8.8 × 1013
10
0.50
17.6
27.1
0.57

Helium400
3.3 × 1014
20
0.94
8.1
1.6
1.6

Table 4.1: Comparison of parameters for hydrogen driven by near-IR and mid-IR lasers and forhelium driven by a 400 nm pulse. A single example intensity is chosen from the range of eachset. We note that the similarity of dimensionless parameters for the 800 and 400 nm cases was bydesign.

152

Figure 4.21: Convergence of macroscopic on-axis radiation with number of cells for (top) Gouy
phase and (bottom) focal phase with g0 = −2, generated by a sin2 pulse of 20 optical cycles
at wavelength of 800 nm and peak intensity of 7.0 × 1013 W/cm2 . Vertical dashed lines denotethe ionization potential (left line) and HHG cut-off (right line). Note different cell numbers thanFig. 4.10.

153
4.6.1

Hydrogen at 1600 nm

HHG driven by mid-infrared pulses has received growing interest in recent years [103–105].
Due to the scaling of the HHG cut-off (IP +3.17 4ωI 2 ) and lower probabilities of ionization (i.e., higher
saturation intensities), longer driving wavelengths allow for production of bright high-harmonicsupercontinua extending into keV energies [31]. Much of the research on mid-IR HHG focuses onhigh driving laser intensities and the resulting ultrahigh harmonics. Here, we seek to investigate
the near-threshold response, where low photon energies and large UP mean that many states crossthrough resonance. As previously mentioned, this spectral region cannot be approximated withthe same approximation methods often used for plateau harmonics. Meanwhile, TDSE calculations
become significantly more computationally expensive due to scaling of pulse duration (∝ λ) andquiver radius (∝ λ2 ) as well as generally higher convergence requirements due to the relatively low
intensity of the HHG signal relative to the fundamental. Efficient interpolation and macroscopic
summation is thus of especial importance at longer wavelengths.
As mentioned above, we have found that longer wavelengths require finer intensity samplingto achieve commensurate convergence, and furthermore that macroscopic summation requires morecells. Therefore, the results of this Section comprise 1, 000 microscopic TDSE simulations at intensity sampling of dI = 1011 W/cm2 , and 107 macroscopic cells. We note that improvements to themethod will be needed in order to expand 1600 nm investigations to inclusion of focal phase – thecompounded convergence requirements exceed current limits of feasibility.
In Fig. 4.22, we present a comparison of a microscopic spectrum (dot-dashed black line) andcalculations where the interpolated single-atom TDSE results are either added up coherently (solidP
red line, as Eq. (4.8)) or incoherently (dashed blue line, with P (ω) = j |aj (ω)|2 ). As for 800
nm, the lack of coherence of the radiation strongly suppresses radiation in the macroscopic result.
Unlike for the lower wavelength, the overall relative intensity of the plateau is greatly reduced, by
about six orders of magnitude compared to the microscopic signal. Additionally, the suppression
does not clearly select odd-harmonic radiation ((2n + 1)ω, denoted by vertical gridlines) compared154
to peaks at other energies, except for radiation above the effective cut-off near 85ω (compared tomicroscopic prediction of 104ω). It thus appears that coherence is poor for all radiation in theplateau, compared to the on- vs. off-harmonic coherence differences seen at lower wavelengths.
This is further supported by Fig. 4.23, where spectra resolved by angle (top) and scanning peakintensity (bottom) show similar lack of harmonic selectivity.

4.6.2

Helium at 400 nm

While atomic hydrogen serves as a straightforward and common target in theory, it is important to also consider how investigations may be extended to experimentally-accessible targets.
Helium is a natural next step as both a common and useful HHG target [31] and the next simplest
atom for simulations. In fact, the lack of inner orbitals is beneficial to both purposes – eliminating
re-absorption of high-energy radiation and avoiding extra steps needed to address non-physicaltransitions to lower occupied states (see Sec. 2.2.3). In our calculations, we have used the helium
SAE potential developed in this thesis (see Ch. 2).
Given the higher ionization potential of helium, we chose a lower wavelength (400 nm) andhigher intensities (by a factor of about 8). This results in similar dimensionless parameters as forhydrogen driven by 800 nm lasers, as seen in Table 4.1. Interpolation and macroscopic summation
are performed in the same manner as for hydrogen, with the same gas jet and laser shape parameters.
While the 400 nm interpolation appears to be converged at higher intensity spacing (≈ 8 × 1012
W/cm2 , consistent with our rule of thumb), for the macroscopic simulations we use the same spacing
as for hydrogen, i.e., dI = 1012 W/cm2 – the TDSE calculations had already been completed forthe finer spacing and the later steps are not significantly impacted by including more points forinterpolation.
Overall, we observe similar features for helium driven by 400 nm lasers as previously observedin the hydrogen results at 800 nm. In Fig. 4.24, on-axis radiation is suppressed between harmonics
by coherent macroscopic summation, but significant off-harmonic peaks do remain. These features
again show, in Fig. 4.25, wider angular dependence than harmonic radiation.

155

Figure 4.22: Comparison of (dot-dash black) single-atom results versus (dashed blue) incoherentand (solid red) coherent macroscopic summation, for HHG in a gas jet of hydrogen atoms driven
by a 1600 nm laser pulse of 10 optical cycles with peak intensity 8.8 × 1013 W/cm2 . Note that theintensity of the spectra are scaled by the signal at the fundamental frequency.

156

Figure 4.23: Macroscopic radiation spectra for a gas jet of hydrogen atoms driven by a 1600 nm
laser. (top) Angle-resolved spectrum for a peak intensity of 8.8 × 1013 W/cm2 . (bottom) On-axis
spectra shown as a function of peak intensity.

157

Figure 4.24: Comparison of single-atom TDSE results and macroscopic TDSE interpolation resultsfor high harmonic generation in helium atoms. On the left, the (a) single-atom and (b) macroscopic
on-axis spectrum is shown as a function of peak intensity. The spectra on the right show thecomparison (single atom: black dashed line, macroscopic: red solid line) at peak intensities (c)I0 = 3.3 × 1014 W/cm2 and (d) I0 = 8.0 × 1013 W/cm2 . (cf. hydrogen results in Fig. 4.11)

158

Figure 4.25: Angle-resolved near-threshold part of the macroscopic radiation spectrum for a gas
jet of helium atoms driven by a 400 nm laser with peak intensity 3.3 × 1014 W/cm2 . (cf. hydrogen
results in Fig. 4.13)

159
One difference appears in the spectrum near 9ω, Fig. 4.26. The spectrum for helium has somesimilar off-harmonic structures as a function of intensity in microscopic simulations (top panel).
Unlike in hydrogen, however, one of these features remains dominant after macroscopic summation(bottom panel). The temporal profile intensity scans (Fig. 4.27) suggest that the same mechanism
is responsible as discussed for hydrogen in Ch. 3.

4.7

Summary
To summarize, we have introduced an efficient method of generating accurate microscopic

single-atom radiation data at arbitrary intensities, based on the results of a relatively small numberof actual time-dependent Schrödinger equation calculations. We have shown that this interpolation
scheme produces spectra matching ab-initio calculations across the entire high-order harmonic
spectra, including the near-threshold regime as well as the entire plateau up to the cut-off. In
the present applications, the error is found to scale with the fourth order of the intensity spacingof the initial TDSE calculations, as expected from the cubic spline method used. While we havefocused on usage of the method for the hydrogen atom, it can be extended to other targets (e.g.,more complex atoms, molecules, solids) as well as other calculation methods (e.g., time-dependent
density functional theory).
We have then applied the interpolation method to compute the radiation response of macroscopic targets. Using a previously tested method for summing the response of microscopic radiators,
we incorporated the efficiently-generated microscopic ab-initio results. This eliminates the limitations of approximative methods in the low-energy portion of the spectrum, near and below theionization threshold. To demonstrate the reliability in using the interpolation method for accurate
macroscopic HHG simulations, we have used the distribution of atoms on a uniform lattice andin a gas jet as test cases. We have shown that angle-resolved macroscopic radiation spectra canbe achieved based on just one hundred ab-initio TDSE calculations. For the uniform lattice, it isshown that the numerical results match the analytical predictions for diffraction patterns.
We have used the TDSE interpolation method to explore features of near-threshold radiation

160

Figure 4.26: Comparison of (top) single-atom results and (bottom) on-axis macroscopic TDSE
interpolation results for near-threshold HHG spectra as a function of peak intensity for a 400 nm,
20 o.c. pulse interacting with helium. Results for each intensity have been normalized to the
maximum signal in the region between 8.75ω and 9.5ω. Red dashed lines denote field-free excitedstate transition energies, and vertical lines identify select intensities: 2.4 × 1014 , 3.3 × 1014 , and
8.0 × 1014 W/cm2 . (cf. hydrogen results in Fig. 4.15)

161

Figure 4.27: Comparison of (top) single-atom results and (bottom) on-axis macroscopic TDSE
interpolation results for temporal profile of radiation within 1ω of the 9th harmonic as a functionof driving pulse intensity in helium for a 400 nm, 20 o.c. pulse. For each intensity, the results havebeen normalized to a maximum of 1 for sake of comparison. A vertical dashed line denotes thepeak of the laser and horizontal dashed lines identify select intensities: 2.4 × 1014 , 3.3 × 1014 , and
8.0 × 1014 W/cm2 . (cf. hydrogen results in Fig. 4.16)

162
in the macroscopic response. Most of the off-harmonic radiation features, which are present at thesame strength as harmonic emission in typical single-atom calculations, are strongly suppressed in
the coherent build-up of the macroscopic signals. The overall macroscopic spectra closely resemble
those obtained in experiments, showing a strong emission at the odd harmonics of the driving laser
frequency and a few features such as weak shoulders close to the harmonics in the near-thresholdregion of the spectrum (incl. the off-harmonic red- and blue-shifted structures studied in Ch. 3– shown to arise from resonant enhancement via a Stark-shifted excited state). We have further
demonstrated that these weak features have a different, i.e., wider, angular spread than the oddharmonic radiation, in agreement with previous experimental and theoretical studies. The temporaland angular dependence both may be useful tools for tailoring pulses.
Investigating an additional focal phase term for highly-focused lasers, we find it has a large
impact on both the angular structure of harmonic radiation and its intensity. Specifically, we haveidentified interference patterns in the angular distributions of below- and near-threshold harmonics.
The interference effects are found in case of a spatial phase distribution for broadband Gaussian
pulses with a negative Porras factor, while they are not present for the monochromatic Gouy phase
distribution. Our analysis indicates that the interferences are due to off-center contributions thatare in-phase with those at the central points in the focus.
Finally, we have briefly explored radiation in two other parameter regimes: hydrogen driven
at longer-wavelengths of a mid-infrared laser, and helium driven by a shorter wavelength laser. For
helium driven by 400 nm laser light, we see behavior consistent with that observed for hydrogenat 800 nm. Future work should in particular investigate intensities near 2 × 1014 W/cm2 , where
off-harmonic radiation remains stronger than 9ω radiation after macroscopic summation. For the
mid-IR results, efficient methods are crucial for overcoming the increasing computational complexity
and costs. At this wavelength, macroscopic radiation in the HHG plateau is strongly suppressed
relative to the microscopic signal, and does not demonstrate the selectivity to odd harmonics foundat lower wavelengths.
Overall, the results obtained for the near-threshold region are in agreement with previous

163
results based on more sophisticated calculations, showing the reliability of the present method for
efficiently exploring a parameter regime rich with complex features.

Chapter 5
Summary and Perspectives

This thesis investigated near-threshold high-order harmonic generation (HHG) through avariety of theoretical techniques. In Chapter 1, we provided background including a simple modelof HHG ignoring the effects of target characteristics, numerical techniques for ab-initio simulationof the time-dependent Schrödinger equation (TDSE) for a single electron in an atomic potential
and the electric field of a laser pulse, and strategies to analyze radiation data.
Next, in Chapter 2 we presented single-active-electron (SAE) potentials for atoms, ions, andmolecules. We demonstrated the utility of these potentials for HHG as well as nonlinear polarizationand ionization studies. Good agreement has been shown between results from separate SAE-TDSE
and time-dependent density functional theory (TDDFT) simulations for the full HHG spectrum,
driven by intense lasers across a variety of wavelengths. This comparison serves to dispute along-standing misconception in the field, and provide a foundation for investigation of impactsof electron correlation on HHG. Additionally, we outlined a strategy to resolve a limitation of
SAE-TDSE simulations in atoms where multiple orbitals may respond to the driving laser.
In Chapter 3, we investigated resonance enhancement of near-threshold HHG. We identified
similar resonance effects in broad parameter regimes, including hydrogen driven by near- and midinfrared pulses, as well as helium with 400 nm lasers. The behavior that may be explicable throughsemi-classical trajectory models (generally at higher intensity) has been differentiated from featuresthat are not consistent with these models. We interpreted results through field-shifted resonanceconditions and dimensionless parameters, outlining two possible regimes where the effect breaks

165
down: first, when a large number of resonances contribute, resulting in noisier signals and second,
when perturbative and/or strongly multiphoton behavior becomes more dominant (small UP /ω
and UP /IP , respectively).
Finally, in Chapter 4 we introduced a method to combine efficient macroscopic simulation
with ab-initio atomic methods. Specifically, we demonstrated that the intensity and CEP dependence of HHG across a wide range of parameters can be efficiently and accurately approximatedthrough simple interpolation. Macroscopic summation then proceeds via the Discrete Dipole Approximation – a method with a proven track record in simulation of higher harmonics within thestrong-field approximation – requiring only hundreds of numerical TDSE simulations. With thisTDSE-DDA method, we showed that near the ionization threshold, off-harmonic radiation featuresare suppressed relative to harmonic radiation, but do remain at significant intensity and with awider angular spread. Using the SAE potentials developed in this thesis, we demonstrated similar
structures in helium driven by a 400 nm laser. The TDSE-DDA method can also be extended to
mid-infrared lasers, where initial investigation showed dramatic suppression of HHG across driving
intensities due to poor coherence – at least at the chosen set of laser and gas parameters. Finally,
we considered an additional spatial phase term present in highly-focused lasers, resulting in a clear
interference structure in the angular dependence of HHG.
Each of these investigations may serve as foundation for future work on understanding HHG.
Therefore, we placed an emphasis throughout on explanations of process, obstacles to overcome,and future directions. For SAE potentials, there are two clear directions for the future: further
careful comparison of SAE-TDSE and TDDFT results to identify dynamic multielectron effects;
and, development and use of additional molecular potentials. With resonance HHG, we have identified similar effects in broad parameter regimes, but placed the majority of our focus on a feature– the impact of the hydrogen 9ω → 3p resonance on radiation near the ninth harmonic – withthe distinction of producing radiation at multiple delays from the resonance. Further investigation
could seek similar cases to illuminate the mechanism, further investigate the parameter boundaries
between forms of HHG (trajectory-like resonance-enhanced, other resonance-enhanced, “normal”),

166
or turn to applications such as pulse-shaping or spectroscopy. With the TDSE-DDA method wehave presented an efficient tool which may easily be extended to a large variety of investigations.
These include observations of interesting effects in HHG from mid-IR pulses and those with additional focal phase terms, which certainly require further attention. Variations of target parameterssuch as gas composition, lattice construction, or location relative to the laser focus may reveal novel
effects in the near-threshold regime. Extensions to more exotic single-radiator simulations (e.g.,molecular targets or bi-elliptical driving fields) hold a great deal of promise. Finally, we have only
just begun to understand the combinations of these three topics, with the helium resonance effects
already showing particular promise in macroscopic simulations.

Bibliography[1] 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).
[2] A. Kramida, Y. Ralchenko, J. Reader, and NIST ASD Team, NIST atomic spectra database(version 5.6.1) (2018).
[3] J. Venzke, R. Reiff, Z. Xue, A. Jaroń-Becker, and A. Becker, Angular momentum distribution
in Rydberg states excited by a strong laser pulse, Phys. Rev. A 98, 043434 (2018).
[4] T. Popmintchev, M.-C. Chen, P. Arpin, M. M. Murnane, and H. C. Kapteyn, The attosecond
nonlinear optics of bright coherent x-ray generation, Nature Photonics 4, 822–832 (2010).
[5] J. H. Heinz, R. Reiff, T. Joyce, A. Jaroń-Becker, A. Becker, and M. Kolesik, Nonlinear
polarization and ionization in O2 : metastable electronic state model, Opt. Express 28, 25022–
25036 (2020).
[6] R. Reiff, J. Venzke, A. Jaroń-Becker, and A. Becker, Interference effects in high harmonic
generation induced by focal phase distribution, submitted to Optics Letters (2021).
[7] H. Hertz, Ueber einen Einfluss des ultravioletten Lichtes auf die electrische Entladung, Annalen der Physik 267, 983–1000 (1887).
[8] A. Einstein, Über einen die Erzeugung und Verwandlung des Lichtes betreffenden heuristischen Gesichtspunkt, Annalen der Physik 4, 132–148 (1905).
[9] M. Göppert-Mayer, Über Elementarakte mit zwei Quantensprüngen, Annalen der Physik
401, 273–294 (1931).
[10] W. Kaiser and C. G. B. Garrett, Two-photon excitation in CaF2 : Eu2+ , Phys. Rev. Lett. 7,
229–231 (1961).
[11] L. Young, K. Ueda, M. Gühr, P. H. Bucksbaum, M. Simon, S. Mukamel, N. Rohringer, K. C.
Prince, C. Masciovecchio, M. Meyer, A. Rudenko, D. Rolles, C. Bostedt, M. Fuchs, D. A.
Reis, R. Santra, H. Kapteyn, M. Murnane, H. Ibrahim, F. Légaré, M. Vrakking, M. Isinger,
D. Kroon, M. Gisselbrecht, A. L’Huillier, H. J. Wörner, and S. R. Leone, Roadmap of ultrafast
x-ray atomic and molecular physics, J. Phys. B: At. Mol. Opt. Phys. 51, 032003 (2018).

168[12] K. Ueda, E. Sokell, S. Schippers, F. Aumayr, H. Sadeghpour, J. Burgdörfer, C. Lemell, X.-M.
Tong, T. Pfeifer, F. Calegari, A. Palacios, F. Martin, P. Corkum, G. Sansone, E. V. Gryzlova,
A. N. Grum-Grzhimailo, M. N. Piancastelli, P. M. Weber, T. Steinle, K. Amini, J. Biegert,
N. Berrah, E. Kukk, R. Santra, A. Müller, D. Dowek, R. R. Lucchese, C. W. McCurdy,
P. Bolognesi, L. Avaldi, T. Jahnke, M. S. Schöffler, R. Dörner, Y. Mairesse, L. Nahon,
O. Smirnova, T. Schlathölter, E. E. B. Campbell, J.-M. Rost, M. Meyer, and K. A. Tanaka,
Roadmap on photonic, electronic and atomic collision physics: I. light–matter interaction, J.
Phys. B: At. Mol. Opt. Phys. 52, 171001 (2019).
[13] A. McPherson, G. Gibson, H. Jara, U. Johann, T. S. Luk, I. A. McIntyre, K. Boyer, andC. K. Rhodes, Studies of multiphoton production of vacuum-ultraviolet radiation in the raregases, J. Opt. Soc. Am. B 4, 595 (1987).
[14] A. L’Huillier, P. Balcou, S. Candel, K. J. Schafer, and K. C. Kulander, Calculations ofhigh-order harmonic-generation processes in xenon at 1064 nm, Phys. Rev. A 46, 2778–2790(1992).
[15] J. Xu, C. I. Blaga, P. Agostini, and L. F. DiMauro, Time-resolved molecular imaging, J.
Phys. B: At. Mol. Opt. Phys. 49, 112001 (2016).
[16] A. S. Alnaser and I. V. Litvinyuk, Subfemtosecond directional control of chemical processes
in molecules, J. Phys. B: At. Mol. Opt. Phys. 50, 032002 (2017).
[17] J. P. Marangos, Development of high harmonic generation spectroscopy of organic molecules
and biomolecules, J. Phys. B: At. Mol. Opt. Phys. 49, 132001 (2016).
[18] H. Yun, S. J. Yun, G. H. Lee, and C. H. Nam, High-harmonic spectroscopy of aligned
molecules, J. Phys. B: At. Mol. Opt. Phys. 50, 022001 (2017).
[19] P. Peng, C. Marceau, and D. M. Villeneuve, Attosecond imaging of molecules using highharmonic spectroscopy, Nature Reviews Physics 1, 144–155 (2019).
[20] H. Reiss, Theoretical methods in quantum optics: S-matrix and Keldysh techniques forstrong-field processes, Prog. Quantum. Electron. 16, 1 – 71 (1992).
[21] M. V. Ammosov, N. B. Delone, and V. P. Krainov, Tunnel ionization of complex atoms andatomic ions in electromagnetic field, in High intensity laser processes, Vol. 664 (International
Society for Optics and Photonics, 1986) pp. 138–141.
[22] 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).
[23] W. Becker, S. P. Goreslavski, D. B. Milošević, and G. G. Paulus, The plateau in abovethreshold ionization: the keystone of rescattering physics, J. Phys. B: At. Mol. Opt. Phys.
51, 162002 (2018).
[24] D. B. Milošević, G. G. Paulus, D. Bauer, and W. Becker, Above-threshold ionization by
few-cycle pulses, J. Phys. B: At. Mol. Opt. Phys. 39, R203–R262 (2006).
[25] N. B. Delone and V. P. Krainov, AC Stark shift of atomic energy levels, Physics-Uspekhi 42,
669–687 (1999).

169[26] 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–2211 (1989).
[27] F. He, C. Ruiz, A. Becker, and U. Thumm, Attosecond probing of instantaneous AC Stark
shifts in helium atoms, J. Phys. B: At. Mol. Opt. Phys. 44, 211001 (2011).
[28] B. Piraux, F. Mota-Furtado, P. F. O’Mahony, A. Galstyan, and Y. V. Popov, Excitation of
Rydberg wave packets in the tunneling regime, Phys. Rev. A 96, 043403 (2017).
[29] 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).
[30] 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: At. Mol. Opt. Phys.
47, 204019 (2014).
[31] T. Popmintchev, M.-C. Chen, D. Popmintchev, P. Arpin, S. Brown, S. Ališauskas, G. Andriukaitis, T. Balčiunas, O. D. Mücke, A. Pugzlys, et al., Bright coherent ultrahigh harmonics
in the keV x-ray regime from mid-infrared femtosecond lasers, Science 336, 1287–1291 (2012).
[32] M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’huillier, and P. B. Corkum, Theory of highharmonic generation by low-frequency laser fields, Phys. Rev. A 49, 2117 (1994).
[33] P. B. Corkum, Plasma perspective on strong field multiphoton ionization, Phys. Rev. Lett.
71, 1994–1997 (1993).
[34] K. Kulander, K. Schafer, and J. Krause, Super-intense laser-atom physics, NATO Advanced
Science Institutes Series 316, 95 (1993).
[35] O. E. Alon, V. Averbukh, and N. Moiseyev, Selection rules for the high harmonic generation
spectra, Phys. Rev. Lett. 80, 3743–3746 (1998).
[36] S. H. Glenzer, L. B. Fletcher, E. Galtier, B. Nagler, R. Alonso-Mori, B. Barbrel, S. B. Brown,
D. A. Chapman, Z. Chen, C. B. Curry, F. Fiuza, E. Gamboa, M. Gauthier, D. O. Gericke,
A. Gleason, S. Goede, E. Granados, P. Heimann, J. Kim, D. Kraus, M. J. MacDonald,
A. J. Mackinnon, R. Mishra, A. Ravasio, C. Roedel, P. Sperling, W. Schumaker, Y. Y. Tsui,
J. Vorberger, U. Zastrau, A. Fry, W. E. White, J. B. Hasting, and H. J. Lee, Matter under
extreme conditions experiments at the Linac coherent light source, J. Phys. B: At. Mol. Opt.
Phys. 49, 092001 (2016).
[37] K. M. Dorney, L. Rego, N. J. Brooks, J. San Román, C.-T. Liao, J. L. Ellis, D. Zusin,
C. Gentry, Q. L. Nguyen, J. M. Shaw, et al., Controlling the polarization and vortex chargeof attosecond high-harmonic beams via simultaneous spin–orbit momentum conservation,
Nature Photonics 13, 123–130 (2019).
[38] P.-C. Huang, C. Hernández-Garcı́a, J.-T. Huang, P.-Y. Huang, C.-H. Lu, L. Rego, D. D.
Hickstein, J. L. Ellis, A. Jaron-Becker, A. Becker, et al., Polarization control of isolated
high-harmonic pulses, Nature Photonics 12, 349–354 (2018).

170[39] S. Beaulieu, S. Camp, D. Descamps, A. Comby, V. Wanie, S. Petit, F. Légaré, K. J. Schafer,
M. B. Gaarde, F. Catoire, and Y. Mairesse, Role of excited states in high-order harmonic
generation, Phys. Rev. Lett. 117, 203001 (2016).
[40] L. Verlet, Computer “experiments” on classical fluids. I. thermodynamical properties ofLennard-Jones molecules, Phys. Rev. 159, 98–103 (1967).
[41] J. A. Hostetter, J. L. Tate, K. J. Schafer, and M. B. Gaarde, Semiclassical approaches tobelow-threshold harmonics, Phys. Rev. A 82, 023401 (2010).
[42] E. Shani, Analysis and numerical performance of methods of solving the time independent
Schrödinger equation for simulation in strong-field physics, Master’s thesis, University ofColorado Boulder (2017).
[43] D. Sorensen, Implicitly restarted Arnoldi/Lanczos methods for large scale eigenvalue calculations, in Parallel Numerical Algorithms, edited by D. Keyes, A. Sameh, and V. Venkatakrishnan (Springer, 1997) Chap. 5, pp. 119–165.
[44] O. A. Sharafeddin, D. J. Kouri, and D. K. Hoffman, Time-dependent treatment of scattering: potential referenced and kinetic energy referenced modified Cayley approaches toatom–diatom collisions and time-scale separations, Canadian Journal of Chemistry 70, 686–
692 (1992).
[45] J. Crank and P. Nicolson, A practical method for numerical evaluation of solutions of partial
differential equations of the heat-conduction type, Mathematical Proceedings of the Cambridge Philosophical Society 43, 50–67 (1947).
[46] Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving
nonsymmetric linear systems, SIAM Journal on Scientific and Statistical Computing 7, 856–
869 (1986).
[47] C. A. Nicolaides and D. R. Beck, The variational calculation of energies and widths of resonances, Phys. Lett. A 65, 11–12 (1978).
[48] B. Simon, The definition of molecular resonance curves by the method of exterior complex
scaling, Phys. Lett. A 71, 211–214 (1979).
[49] 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).
[50] L. Keldysh, Ionization in the field of a strong electromagnetic wave, JETP 20, 1307 (1965).
[51] F. H. M. Faisal, Multiple absorption of laser photons by atoms, J. Phys. B: At. Mol. Opt.
Phys. 6, L89–L92 (1973).
[52] H. R. Reiss, Effect of an intense electromagnetic field on a weakly bound system, Phys. Rev.
A 22, 1786–1813 (1980).
[53] J. A. Pérez-Hernández, L. Roso, and L. Plaja, Harmonic generation beyond the strong-fieldapproximation: the physics behind the short-wave-infrared scaling laws, Opt. Express 17,9891–9903 (2009).

171[54] K. Amini, J. Biegert, F. Calegari, A. Chacón, M. F. Ciappina, A. Dauphin, D. K. Efimov,
C. F. de Morisson Faria, K. Giergiel, P. Gniewek, A. S. Landsman, M. Lesiuk, M. Mandrysz,
A. S. Maxwell, R. Moszyński, L. Ortmann, J. A. Pérez-Hernández, A. Picón, E. Pisanty,
J. Prauzner-Bechcicki, K. Sacha, N. Suárez, A. Zaı̈r, J. Zakrzewski, and M. Lewenstein,
Symphony on strong field approximation, Reports on Progress in Physics 82, 116001 (2019).
[55] A. V. Oppenheim, R. W. Schafer, and J. R. Buck, Filter design techniques, in
Discrete-Time Signal Processing (Prentice Hall, 1999) p. 468–471.
[56] J. M. Lilly and S. C. Olhede, Higher-order properties of analytic wavelets, IEEE Transactions
on Signal Processing 57, 146–160 (2009).
[57] J. M. Lilly and S. C. Olhede, Generalized Morse wavelets as a superfamily of analytic wavelets,
IEEE Transactions on Signal Processing 60, 6036–6041 (2012).
[58] J. M. Lilly, Element analysis: a wavelet-based method for analysing time-localized events innoisy time series, Proc. R. Soc. A. 473, 20160776 (2017).
[59] A. Crawford-Uranga, U. De Giovannini, E. Räsänen, M. J. T. Oliveira, D. J. Mowbray, G. M.
Nikolopoulos, E. T. Karamatskos, D. Markellos, P. Lambropoulos, S. Kurth, and A. Rubio,
Time-dependent density-functional theory of strong-field ionization of atoms by soft x rays,
Phys. Rev. A 90, 033412 (2014).
[60] X. M. Tong and S.-I. Chu, Time-dependent density-functional theory for strong-field multiphoton processes: Application to the study of the role of dynamical electron correlation in
multiple high-order harmonic generation, Phys. Rev. A 57, 452–461 (1998).
[61] J. J. Carrera and S.-I. Chu, Ab initio time-dependent density-functional-theory study of thefrequency comb structure, coherence, and dephasing of multielectron systems in the VUVXUV regimes via high-order harmonic generation, Phys. Rev. A 79, 063410 (2009).
[62] M. D. Perry, A. Szoke, and K. C. Kulander, Resonantly enhanced above-threshold ionizationof helium, Phys. Rev. Lett. 63, 1058–1061 (1989).
[63] K. C. Kulander and T. N. Rescigno, Effective potentials for time-dependent calculations ofmultiphoton processes in atoms, Comput. Phys. Commun. 63, 523 (1991).
[64] K. C. Kulander, K. J. Schafer, and J. L. Krause, Single-active electron calculation of multiphoton process in krypton, Int. J. Quant. Chem. 40, 415 (1991).
[65] H. G. Muller and F. C. Kooiman, Bunching and focusing of tunneling wave packets in enhancement of high-order above-threshold ionization, Phys. Rev. Lett. 81, 1207 (1998).
[66] H. G. Muller, Numerical simulation of high-order above-threshold-ionization enhancement inargon, Phys. Rev. A 60, 1341 (1999).
[67] M. J. Nandor, M. A. Walker, L. D. Van Woerkom, and H. G. Muller, Detailed comparison ofabove-threshold-ionization spectra from accurate numerical integrations and high-resolution
measurements, Phys. Rev. A 60, R1771 (1999).
[68] X. M. Tong and C. D. Lin, Empirical formula for static field ionization rates of atoms andmolecules by lasers in the barrier-suppression regime, J. Phys. B: At. Mol. Opt. Phys. 38,
2593 (2005).

172[69] M. Awasthi, Y. V. Vanne, A. Saenz, A. Castro, and P. Decleva, Single-active-electron approximation for describing molecules in ultrashort laser pulses and its application to molecularhydrogen, Phys. Rev. A 77, 063403 (2008).
[70] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
[71] W. Kohn and L. J. Sham, One-particle properties of an inhomogeneous interacting electrongas, Phys. Rev. 140, A1133 (1965).
[72] J. B. Krieger, Y. Li, and G. J. Iafrate, Exact relations in the optimized effective potential
method employing an arbitrary Exc [ψiσ ], Phys. Lett. A 148, 470 (1990).
[73] X. M. Tong and S.-I. Chu, Density-functional theory with optimized effective potential andself-interaction correction for ground states and autoionizing resonances, Phys. Rev. A 55,
3406–3416 (1997).
[74] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations
for many-electron systems, Phys. Rev. B 23, 10 (1981).
[75] M. A. L. Marques, A. Castro, G. F. Bertsch, and A. Rubio, octopus: a first-principles toolfor excited electron–ion dynamics, Comput. Phys. Commun. 151, 60 (2003).
[76] A. Castro, H. Appel, M. Oliveira, C. A. Rozzi, X. Andrade, F. Lorenzen, M. A. L. Marques,
E. K. U. Gross, and A. Rubio, Octopus: a tool for the application of time-dependent density
functional theory, Phys. Stat. Sol. B 243, 2465 (2006).
[77] X. Andrade, D. A. Strubbe, U. De Giovannini, A. H. Larsen, M. J. T. Oliveira, J. AlberdiRodriguez, A. Varas, I. Theophilou, N. Helbig, M. Verstraete, L. Stella, F. Nogueira,
A. Aspuru-Guzik, A. Castro, M. A. L. Marques, and A. Rubio, Real-space grids and theOctopus code as tools for the development of new simulation approaches for electronic systems, Phys. Chem. Chem. Phys. 17, 31371 (2015).
[78] M. A. Marques, M. J. Oliveira, and T. Burnus, Libxc: A library of exchange and correlation
functionals for density functional theory, Comput. Phys. Commun. 183, 2272–2281 (2012).
[79] D. D. Johnson, Modified Broyden’s method for accelerating convergence in self-consistentcalculations, Phys. Rev. B 38, 12807–12813 (1988).
[80] R. Kosloff and D. Kosloff, Absorbing boundaries for wave propagation problems, J. Comput.
Phys. 63, 363–376 (1986).
[81] K. Levenberg, A method for the solution of certain nonlinear problems, Q. Appl. Math. 2,
164 (1944).
[82] D. W. Marquardt, An algorithm for least-squares estimation of nonlinear parameters, SIAM
J. Appl. Math. 11, 431 (1963).
[83] M. Newville, A. Nelson, T. Stensitzki, A. Ingargiola, and D. Allan, LMFIT: Non-linear leastsquare minimization and curve-fitting for python (2014).
[84] A. D. Bandrauk, S. Chelkowski, D. J. Diestler, J. Manz, and K.-J. Yuan, Quantum simulationof high-order harmonic spectra of the hydrogen atom, Phys. Rev. A 79, 023403 (2009).

173[85] 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).
[86] A. D. Shiner, B. E. Schmidt, C. Trallero-Herrero, H. J. Wörner, S. Patchkovskii, P. B.
Corkum, J. C. Kieffer, F. Legare, and D. M. Villeneuve, Probing collective multi-electron
dynamics in xenon with high-harmonic spectroscopy, Nat. Phys. 7, 2011 (2011).
[87] S. Pabst, L. Greenman, D. A. Mazziotti, and R. Santra, Impact of multichannel and multipole
effects on the Cooper minimum in the high-order-harmonic spectrum of argon, Phys. Rev. A
85, 023411 (2012).
[88] D. Facciala, S. Pabst, B. D. Bruner, A. G. Ciriolo, S. De Silvestri, M. Devetta, M. Negro,
H. Soifer, S. Stagira, N. Dudovich, and C. Vozzi, Probe of multielectron dynamics in xenon
by caustics in high-order harmonic generation, Phys. Rev. Lett. 117, 093902 (2016).
[89] A. W. Bray, D. Freeman, S. Eckart, and A. S. Kheifets, Correlation enhancement of high-orderharmonic generation in Xe, Phys. Rev. A 100, 013404 (2019).
[90] M. Kolesik, J. M. Brown, A. Teleki, P. Jakobsen, J. V. Moloney, and E. M. Wright, Metastableelectronic states and nonlinear response for high-intensity optical pulses, Optica 1, 323–331(2014).
[91] A. Bahl, V. P. Majety, A. Scrinzi, and M. Kolesik, Nonlinear optical response in molecular
nitrogen: from ab-initio calculations to optical pulse simulations, Opt. Lett. 42, 2295–2298(2017).
[92] J. Venzke, A. Jaroń-Becker, and A. Becker, Ionization of helium by an ultrashort extremeultraviolet laser pulse, J. Phys. B: At. Mol. Opt. Phys. 53, 085602 (2020).
[93] L. Kolanz, Anomalous ionization of atoms in the strong-field regime, Undergraduate honors
thesis, Department of Physics, University of Colorado Boulder (2020).
[94] C. F. Pauletti, E. Coccia, and E. Luppi, Role of exchange and correlation in high-harmonic
generation spectra of H2 , N2 , and CO2 : Real-time time-dependent electronic-structure approaches, J. Chem. Phys. 154, 014101 (2021).
[95] D. C. Yost, T. R. Schibli, J. Ye, J. L. Tate, J. Hostetter, M. B. Gaarde, and K. J. Schafer,
Vacuum-ultraviolet frequency combs from below-threshold harmonics, Nature Physics 5, 815–
820 (2009).
[96] M. Chini, X. Wang, Y. Cheng, H. Wang, Y. Wu, E. Cunningham, P.-C. Li, J. Heslar, D. A.
Telnov, S.-I. Chu, et al., Coherent phase-matched VUV generation by field-controlled boundstates, Nature Photonics 8, 437–441 (2014).
[97] 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).
[98] S. Bengtsson, E. W. Larsen, D. Kroon, S. Camp, M. Miranda, C. Arnold, A. L’Huillier,K. Schafer, M. Gaarde, L. Rippe, and J. Mauritsson, Space–time control of free induction
decay in the extreme ultraviolet, Nature Photonics 11, 252–258 (2017).

174[99] W.-H. Xiong, J.-Z. Jin, L.-Y. Peng, and Q. Gong, Numerical observation of two sets of loworder harmonics near the ionization threshold, Phys. Rev. A 96, 023418 (2017).
[100] S. Camp, S. Beaulieu, K. J. Schafer, and M. B. Gaarde, Resonantly-initiated quantum trajectories and their role in the generation of near-threshold harmonics, J. Phys. B: At. Mol.
Opt. Phys. 51, 064001 (2018).
[101] H. Yun, J. H. Mun, S. I. Hwang, S. B. Park, I. A. Ivanov, C. H. Nam, and K. T. Kim,
Coherent extreme-ultraviolet emission generated through frustrated tunnelling ionization,
Nature Photonics 12, 620–624 (2018).
[102] M. R. Miller, C. Hernández-Garcı́a, A. Jaroń-Becker, and A. Becker, Targeting multiple
rescatterings through VUV-controlled high-order-harmonic generation, Phys. Rev. A 90,
053409 (2014).
[103] C. Hernández-Garcı́a, C. G. Durfee, D. D. Hickstein, T. Popmintchev, A. Meier, M. M.
Murnane, H. C. Kapteyn, I. J. Sola, A. Jaron-Becker, and A. Becker, Schemes for generationof isolated attosecond pulses of pure circular polarization, Phys. Rev. A 93, 043855 (2016).
[104] K. Werner, M. G. Hastings, A. Schweinsberg, B. L. Wilmer, D. Austin, C. M. Wolfe, M. Kolesik, T. R. Ensley, L. Vanderhoef, A. Valenzuela, and E. Chowdhury, Ultrafast mid-infrared
high harmonic and supercontinuum generation with n2 characterization in zinc selenide, Opt.
Express 27, 2867–2885 (2019).
[105] C. Hernández-Garcı́a, T. Popmintchev, M. M. Murnane, H. C. Kapteyn, L. Plaja, A. Becker,
and A. Jaron-Becker, Isolated broadband attosecond pulse generation with near- and midinfrared driver pulses via time-gated phase matching, Opt. Express 25, 11855–11866 (2017).
[106] C. Hernández-Garcı́a, J. A. Pérez-Hernández, J. Ramos, E. C. Jarque, L. Roso, and L. Plaja,
High-order harmonic propagation in gases within the discrete dipole approximation, Phys.
Rev. A 82, 033432 (2010).
[107] T. Popmintchev, M.-C. Chen, A. Bahabad, M. Gerrity, P. Sidorenko, O. Cohen, I. P. Christov,M. M. Murnane, and H. C. Kapteyn, Phase matching of high harmonic generation in the softand hard x-ray regions of the spectrum, Proceedings of the National Academy of Sciences
106, 10516–10521 (2009).
[108] E. Lorin, S. Chelkowski, and A. Bandrauk, A numerical Maxwell–Schrödinger model for
intense laser–matter interaction and propagation, Comput. Phys. Commun. 177, 908 – 932(2007).
[109] M. B. Gaarde, J. L. Tate, and K. J. Schafer, Macroscopic aspects of attosecond pulse generation, J. Phys. B: At. Mol. Opt. Phys. 41, 132001 (2008).
[110] E. Lorin, M. Lytova, A. Memarian, and A. Bandrauk, Development of nonperturbative nonlinear optics models including effects of high order nonlinearities and of free electron plasma:
Maxwell–Schrödinger equations coupled with evolution equations for polarization effects, andthe SFA-like nonlinear optics model, J. Phys. A Math. Theor. 48, 105201 (2015).

175[111] V. V. Strelkov, V. T. Platonenko, A. F. Sterzhantov, and M. Y. Ryabikin, Attosecond electromagnetic pulses: generation, measurement, and application. Generation of high-order harmonics of an intense laser field for attosecond pulse production, Physics-Uspekhi 59, 425–445(2016).
[112] I. P. Christov, M. M. Murnane, and H. C. Kapteyn, Generation and propagation of attosecond
x-ray pulses in gaseous media, Phys. Rev. A 57, R2285–R2288 (1998).
[113] M. Geissler, G. Tempea, A. Scrinzi, M. Schnürer, F. Krausz, and T. Brabec, Light propagationin field-ionizing media: Extreme nonlinear optics, Phys. Rev. Lett. 83, 2930–2933 (1999).
[114] E. Priori, G. Cerullo, M. Nisoli, S. Stagira, S. De Silvestri, P. Villoresi, L. Poletto, P. Ceccherini, C. Altucci, R. Bruzzese, and C. de Lisio, Nonadiabatic three-dimensional model ofhigh-order harmonic generation in the few-optical-cycle regime, Phys. Rev. A 61, 063801(2000).
[115] N. H. Shon, A. Suda, and K. Midorikawa, Generation and propagation of high-order harmonics
in high-pressure gases, Phys. Rev. A 62, 023801 (2000).
[116] M. B. Gaarde and K. J. Schafer, Space-time considerations in the phase locking of highharmonics, Phys. Rev. Lett. 89, 213901 (2002).
[117] V. Tosa, H. T. Kim, I. J. Kim, and C. H. Nam, High-order harmonic generation by chirped
and self-guided femtosecond laser pulses. I. Spatial and spectral analysis, Phys. Rev. A 71,
063807 (2005).
[118] T. Ruchon, C. P. Hauri, K. Varjú, E. Mansten, M. Swoboda, R. López-Martens, and
A. L’Huillier, Macroscopic effects in attosecond pulse generation, New Journal of Physics
10, 025027 (2008).
[119] V. Tosa, K. T. Kim, and C. H. Nam, Macroscopic generation of attosecond-pulse trains instrongly ionized media, Phys. Rev. A 79, 043828 (2009).
[120] C. Jin, A.-T. Le, and C. D. Lin, Retrieval of target photorecombination cross sections fromhigh-order harmonics generated in a macroscopic medium, Phys. Rev. A 79, 053413 (2009).
[121] E. Lorin, S. Chelkowski, E. Zaoui, and A. Bandrauk, Maxwell–Schrödinger–plasma (MASP)
model for laser–molecule interactions: Towards an understanding of filamentation with intense ultrashort pulses, Physica D: Nonlinear Phenomena 241, 1059 – 1071 (2012).
[122] E. Lorin and A. Bandrauk, Efficient and accurate numerical modeling of a micro–macro
nonlinear optics model for intense and short laser pulses, Journal of Computational Science3, 159 – 168 (2012).
[123] M. C. Kohler and K. Z. Hatsagortsyan, Macroscopic aspects of relativistic x-ray-assisted
high-order-harmonic generation, Phys. Rev. A 85, 023819 (2012).
[124] M. Kolesik, E. M. Wright, J. Andreasen, J. M. Brown, D. R. Carlson, and R. J. Jones,
Space-time resolved simulation of femtosecond nonlinear light-matter interactions using aholistic quantum atomic model : Application to near-threshold harmonics, Opt. Express 20,
16113–16128 (2012).

176[125] P.-C. Li and S.-I. Chu, Effects of macroscopic propagation on spectra of broadband supercontinuum harmonics and isolated-attosecond-pulse generation: Coherent control of the electron
quantum trajectories in two-color laser fields, Phys. Rev. A 86, 013411 (2012).
[126] S. A. Berman, J. Dubois, C. Chandre, M. Perin, and T. Uzer, Coherent buildup of high-orderharmonic radiation: The classical perspective, Phys. Rev. A 97, 061402 (2018).
[127] C. Hernández-Garcı́a, A. Picón, J. San Román, and L. Plaja, Attosecond extreme ultraviolet
vortices from high-order harmonic generation, Phys. Rev. Lett. 111, 083602 (2013).
[128] D. D. Hickstein, F. J. Dollar, P. Grychtol, J. L. Ellis, R. Knut, C. Hernández-Garcı́a, D. Zusin,
C. Gentry, J. M. Shaw, T. Fan, et al., Non-collinear generation of angularly isolated circularly
polarized high harmonics, Nature Photonics 9, 743–750 (2015).
[129] W.-H. Xiong, L.-Y. Peng, and Q. Gong, Recent progress of below-threshold harmonic generation, J. Phys. B: At. Mol. Opt. Phys. 50, 032001 (2017).
[130] J. B. Watson, A. Sanpera, X. Chen, and K. Burnett, Harmonic generation from a coherent
superposition of states, Phys. Rev. A 53, R1962–R1965 (1996).
[131] E. Toma, P. Antoine, A. de Bohan, and H. Muller, Resonance-enhanced high-harmonic generation, J. Phys. B: At. Mol. Opt. Phys. 32, 5843 (1999).
[132] M. B. Gaarde and K. J. Schafer, Enhancement of many high-order harmonics via a single
multiphoton resonance, Phys. Rev. A 64, 013820 (2001).
[133] R. Taı̈eb, V. Véniard, J. Wassaf, and A. Maquet, Roles of resonances and recollisions instrong-field atomic phenomena. II. High-order harmonic generation, Phys. Rev. A 68, 033403(2003).
[134] D. Milošević, High-energy stimulated emission from plasma ablation pumped by resonant
high-order harmonic generation, J. Phys. B: At. Mol. Opt. Phys. 40, 3367 (2007).
[135] V. Strelkov, Role of autoionizing state in resonant high-order harmonic generation and attosecond pulse production, Phys. Rev. Lett. 104, 123901 (2010).
[136] A. Baltuška, T. Udem, M. Uiberacker, M. Hentschel, E. Goulielmakis, C. Gohle,
R. Holzwarth, V. S. Yakovlev, A. Scrinzi, T. W. Hänsch, and F. Krausz, Attosecond controlof electronic processes by intense light fields, Nature 421, 611–615 (2003).
[137] I. P. Christov, M. M. Murnane, and H. C. Kapteyn, High-Harmonic Generation of Attosecond
Pulses in the “Single-Cycle” Regime, Phys. Rev. Lett. 78, 1251–1254 (1997).
[138] G. G. Paulus, F. Lindner, H. Walther, A. Baltuška, E. Goulielmakis, M. Lezius, and F. Krausz,
Measurement of the Phase of Few-Cycle Laser Pulses, Phys. Rev. Lett. 91, 253004 (2003).
[139] M. F. Kling, C. Siedschlag, A. J. Verhoef, J. I. Khan, M. Schultze, T. Uphues, Y. Ni, M. Uiberacker, M. Drescher, F. Krausz, and M. J. J. Vrakking, Control of Electron Localization in
Molecular Dissociation, Science 312, 246–248 (2006).
[140] M. Krüger, M. Schenk, and P. Hommelhoff, Attosecond control of electrons emitted from ananoscale metal tip, Nature 475, 78–81 (2011).

177[141] D. Guénot, D. Gustas, A. Vernier, B. Beaurepaire, F. Böhle, M. Bocoum, M. Lozano, A. Jullien, R. Lopez-Martens, A. Lifschitz, and J. Faure, Relativistic electron beams driven by kHz
single-cycle light pulses, Nature Photonics 11, 293–296 (2017).
[142] M. A. Porras, Characterization of the electric field of focused pulsed Gaussian beams forphase-sensitive interactions with matter, Optics Letters 34, 1546–1548 (2009).
[143] M. A. Porras, Z. L. Horváth, and B. Major, Three-dimensional carrier-envelope-phase mapof focused few-cycle pulsed Gaussian beams, Phys. Rev. A 98, 063819 (2018).
[144] L. G. Gouy, Sur une propriété nouvelle des ondes lumineuses, Compt. Rendue Acad. Sci. 110,
1251–1253 (1890).
[145] M. F. Parsons, F. Huber, A. Mazurenko, C. S. Chiu, W. Setiawan, K. Wooley-Brown, S. Blatt,and M. Greiner, Site-resolved imaging of Fermionic 6 Li in an optical lattice, Phys. Rev. Lett.
114, 213002 (2015).
[146] S. L. Campbell, R. B. Hutson, G. E. Marti, A. Goban, N. Darkwah Oppong, R. L. McNally,
L. Sonderhouse, J. M. Robinson, W. Zhang, B. J. Bloom, and J. Ye, A Fermi-degenerate
three-dimensional optical lattice clock, Science 358, 90–94 (2017).
[147] D. Jaksch, Optical lattices, ultracold atoms and quantum information processing, Contemporary Physics 45, 367–381 (2004).
[148] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen, and U. Sen, Ultracold atomic
gases in optical lattices: mimicking condensed matter physics and beyond, Advances inPhysics 56, 243–379 (2007).
[149] T. Esslinger, Fermi-Hubbard physics with atoms in an optical lattice, Annual Review ofCondensed Matter Physics 1, 129–152 (2010).
[150] M. Lewenstein, A. Sanpera, and V. Ahufinger, Ultracold atoms in optical lattices: Simulating
quantum many-body systems, Oxford University Press (2012).
[151] C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, Science
357, 995–1001 (2017).
[152] P. Wessels, B. Ruff, T. Kroker, A. K. Kazansky, N. M. Kabachnik, K. Sengstock, M. Drescher,and J. Simonet, Absolute strong-field ionization probabilities of ultracold rubidium atoms,
Commun. Phys. 1, 1–9 (2018).
[153] T. Auguste, P. Salières, A. S. Wyatt, A. Monmayrant, I. A. Walmsley, E. Cormier, A. Zaı̈r,
M. Holler, A. Guandalini, F. Schapper, J. Biegert, L. Gallmann, and U. Keller, Theoreticaland experimental analysis of quantum path interferences in high-order harmonic generation,
Phys. Rev. A 80, 033817 (2009).
[154] F. Schapper, M. Holler, T. Auguste, A. Zaı̈r, M. Weger, P. Salières, L. Gallmann, and
U. Keller, Spatial fingerprint of quantum path interferences in high order harmonic generation, Opt. Express 18, 2987–2994 (2010).
[155] M. Kretschmar, C. Hernández-Garcı́a, D. S. Steingrube, L. Plaja, U. Morgner, andM. Kovačev, Spatial contributions of electron trajectories to high-order-harmonic radiation
originating from a semi-infinite gas cell, Phys. Rev. A 88, 013805 (2013).

178[156] J. Jiang, P. Wei, Z. Zeng, J. Miao, Y. Zheng, X. Ge, C. Li, R. Li, and Z. Xu, Spatially andspectrally resolved quantum-path tracing in high-order harmonic generation, Opt. Express
22, 15975–15981 (2014).
[157] F. Catoire, A. Ferré, O. Hort, A. Dubrouil, L. Quintard, D. Descamps, S. Petit, F. Burgy,
E. Mével, Y. Mairesse, and E. Constant, Complex structure of spatially resolved high-orderharmonic spectra, Phys. Rev. A 94, 063401 (2016).
[158] D. Hoff, M. Krüger, L. Maisenbacher, A. M. Sayler, G. G. Paulus, and P. Hommelhoff, Tracing
the phase of focused broadband laser pulses, Nature Physics 13, 947–951 (2017).
[159] Y. Zhang, D. Zille, D. Hoff, P. Wustelt, D. Würzler, M. Möller, A. M. Sayler, and G. G.
Paulus, Observing the Importance of the Phase-Volume Effect for Few-Cycle Light-Matter
Interactions, Phys. Rev. Lett. 124, 133202 (2020).
[160] S. Beaulieu, S. Camp, D. Descamps, A. Comby, V. Wanie, S. Petit, F. Légaré, K. J. Schafer,
M. B. Gaarde, F. Catoire, and Y. Mairesse, Role of Excited States In High-order Harmonic
Generation, Phys. Rev. Lett. 117, 203001 (2016).
[161] F. Ferray, A. L’Huillier, X. F. Li, L. A. Lompre, G. Mainfray, and C. Manus, Multipleharmonic conversion of 1064 nm radiation in rare gases, J. Phys. B: At. Mol. Opt. Phys. 21,
L31 (1988).
[162] K. C. Kulander, Multiphoton ionization of hydrogen: A time-dependent theory, Phys. Rev.
A 35, 445 (1987).
[163] K. J. Schafer and K. C. Kulander, Energy analysis of time-dependent wave functions: Application to above-threshold ionization, Phys. Rev. A 42, 5794 (1990).
[164] J. S. Parker, B. S. J. Doherty, K. J. Meharg, and K. T. Taylor, Time delay between singly
and doubly ionizing wavepackets in laser-driven helium, J. Phys. B: At. Mol. Opt. Phys. 36,
L393 (2003).
[165] E. Cormier and P. Lambropoulos, Above-threshold ionization spectrum of hydrogen using
B-spline functions, J. Phys. B: At. Mol. Opt. Phys. 30, 77 (1977).
[166] H. Bachau, E. Cormier, P. Decleva, J. E. Hansen, and F. Martin, Applications of B-splines
in atomic and molecular physics, Rep. Prog. Phys. 64, 1815 (2001).
[167] K. C. Kulander, Time-dependent Hartree-Fock theory of multiphoton ionization: Helium,
Phys. Rev. A 36, 2726 (1987).
[168] X. Li, S. M. Smith, A. N. Markevitch, D. A. Romanov, R. J. Levis, and H. B. Schlegel, A timedependent Hartree–Fock approach for studying the electronic optical response of molecules
in intense fields, Phys. Chem. Chem. Phys. 7, 233 (2005).
[169] J. Caillat, J. Zanghellini, M. Kitzler, O. Koch, W. Kreuzer, and A. Scrinzi, Correlated multielectron systems in strong laser fields: A multiconfiguration time-dependent Hartree-Fockapproach, Phys. Rev. A 71, 012712 (2005).
[170] M. Nest, R. Padmanaban, and P. Saalfrank, Time-dependent approach to electronically excited states of molecules with the multiconfiguration time-dependent Hartree-Fock method,
J. Chem. Phys. 126, 214106 (2007).

179[171] S. Sukiasyan, C. McDonald, C. Destefani, M. Y. Ivanov, and T. Brabec, Multielectron correlation in high-harmonic generation: A 2D model analysis, Phys. Rev. Lett. 102, 223002(2009).
[172] D. Haxton and C. McCurdy, Ultrafast population transfer to excited valence levels of amolecule driven by x-ray pulses, Phys. Rev. A 90, 053426 (2014).
[173] M. A. Lysaght, H. W. van der Hart, and P. B. Burke, Time-dependent R-matrix theory for
ultrafast atomic processes, Phys. Rev. A 79, 053411 (2009).
[174] O. Hassouneh, A. C. Brown, and H. W. van der Hart, Harmonic generation by noble-gas
atoms in the near-IR regime using ab initio time-dependent R-matrix theory, Phys. Rev. A
90, 043418 (2014).
[175] P. Krause, T. Klamroth, and P. Saalfrank, Time-dependent configuration-interaction calculations of laser-pulse-driven many-electron dynamics: Controlled dipole switching in lithium
cyanide, J. Chem. Phys. 123, 074105 (2005).
[176] L. Greenman, P. J. Ho, S. Pabst, E. Kamarchik, D. A. Mazziotti, and R. Santra, Implementation of the time-dependent configuration-interaction singles method for atomic strong-fieldprocesses, Phys. Rev. A 82, 023406 (2010).
[177] P. Krause, J. A. Song, and H. B. Schlegel, Strong field ionization rates simulated with timedependent configuration interaction and an absorbing potential, J. Chem. Phys. 140, 174113(2014).
[178] C. Huber and T. Klammroth, Explicitly time-dependent coupled cluster singles doubles calculations of laser-driven many-electron dynamics, J. Chem. Phys. 134, 054113 (2011).
[179] H. Miyagi and L. B. Madsen, Time-dependent restricted-active-space self-consistent-field theory for laser-driven many-electron dynamics, Phys. Rev. A 87, 062511 (2013).
[180] T. Sato and K. L. Ishikawa, Time-dependent complete-active-space self-consistent-field
method for multielectron dynamics in intense laser fields, Phys. Rev. A 88, 023402 (2013).
[181] T. Sato, K. L. Ishikawa, I. Brezinova, F. Lackner, S. Nagele, and J. Burgdörfer, Timedependent complete-active-space self-consistent-field method for atoms: Application to highorder harmonic generation, Phys. Rev. A 94, 023405 (2016).
[182] M. Brics and D. Bauer, Time-dependent renormalized natural orbital theory applied to thetwo-electron spin-singlet case: Ground state, linear response, and autoionization, Phys. Rev.
A 88, 052514 (2013).
[183] M. Brics, J. Rapp, and D. Bauer, Strong-field absorption and emission of radiation in twoelectron systems calculated with time-dependent natural orbitals, Phys. Rev. A 93, 013404(2016)