Perturbative and ab-initio calculations of electricalsusceptibilities of atoms
byAndrew Spott
B.S. Physics, University of Washington, 2010
M.S. Physics, University of Colorado, 2014

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
2017This thesis entitled:
Perturbative and ab-initio calculations of electrical susceptibilities of atoms
written by Andrew Spotthas been approved for the Department of Physics

Prof. Andreas Becker

Prof. Agnieszka Jaroń-Becker

Date

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

Spott, Andrew (Ph.D., Physics)
Perturbative and ab-initio calculations of electrical susceptibilities of atoms
Thesis directed by Prof. Andreas Becker

Perturbative nonlinear optics consists of many powerful predictive theoretical methods, including the perturbative series of observables related to the interaction of light with matter. The
light intensity limits of such series have been studied in the past for highly nonlinear processes
such as above threshold ionization and high harmonic generation. A more recent debate focuses on
the limits of applicability of perturbation theory for the nonlinear electrical susceptibility and thenonlinear index of refraction of atoms, which are important parameters to study, for example, forfilamentation of laser pulses in nonlinear media.
In this thesis we analyze theoretical predictions for the electrical susceptibility of atoms forthe transition from the perturbative to the nonperturbative intensity regime. To this end, we apply
a numerical basis state method that allows us to perform respective calculations in the frameworkof perturbation theory as well as using ab-initio methods. The results let us identify the intensity
at which the application of perturbation theory breaks down. Furthermore, we provide an analysisof the nonlinear susceptibility as a function of time during the interaction with the laser pulse andfind that theoretical predictions are in good agreement with recent experimental data.

To my family: my brother Alex and my parents Tom and Sherry, and to my love Alejandra.

Acknowledgements

No man is an island, and I am no exception. A large part of this thesis has only been
completed due to the support of those around me.
Much of a PhD is learning the process and techniques of research, as well as the research itself.
Towards this, I am in debt to the members of my research group for their help. Dr. Michelle Miller
was an invaluable sounding board of ideas and thoughts, in particular with regards to presenting
and writing the results when they occurred. Cory Goldsmith proved a crucial sounding board
and compatriot when discussing science. Without our discussions, this thesis would likely containmuch more in the way of small misunderstandings and misinterpretations. Towards the end of thisjourney, Brynn Reiff and Joel Venzke have been important partners in the musings on physical lawthat are such an important part of research.
My advisers Andreas Becker and Agnieszka Jaroń-Becker have been instrumental in mysuccess in this endeavor. Their support as scientists – both in the direction of the research andnavigating the unknown in the quest for an answer – is a large part of what made me the scientist
I have become.
I would also like to thank the other members of my thesis committee: Prof. David Jonas,
Prof. Margaret Murnane and Dr. José D’incao. It was a privilege to have them on my committee.
Graduate school is not just about the research however. The mental fortitude to keep returning to solve problems despite setbacks and frustration is made easier with the support of friends
and family. Gregory Petropoulos and Daniel Weingarten deserve special mention in this capacity.
Their friendship turned what might otherwise be a grind into an enjoyable experience, and I amvivery grateful for their friendship.
I would not be here without my parents, who instilled in me the creativity, curiosity, anddrive that got me where I am today. For this I am inordinately grateful.
I would like to thank my brother, Alexander, as well. His skill at questioning basic assumptions has ensured my own understanding of physics.
Finally, these acknowledgements would not be complete without thanking Alejandra. She
has been with me throughout graduate school, and her unfailing optimism and unending support
are supremely appreciated.

Contents

Chapter
1 Introduction1

1.1

Self-focusing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.2

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

5

1.3

Higher order Kerr effect (HOKE) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.4

Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92 Theory

11

2.1

Index of refraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.2

General properties of the electrical susceptibility . . . . . . . . . . . . . . . . . . . . 14

2.3

Perturbation Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.4

2.3.1

Brief review of perturbation theory . . . . . . . . . . . . . . . . . . . . . . . . 16

2.3.2

Dipole moment and single atom response functions . . . . . . . . . . . . . . . 17

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3 Numerical Methods
3.1

3.2

23

Overview of the method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.1

Single active electron models . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.1.2

Separation of variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

Numerical basis state method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.2.1

Representations of differential equations . . . . . . . . . . . . . . . . . . . . . 29

viii

3.3

3.4

3.2.2

Numerical basis state method . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.2.3

Properties of the basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

3.2.4

Implementation of numerical basis state method . . . . . . . . . . . . . . . . 42

Time dependent calculations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

3.3.1

Propagation calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

3.3.2

Absorbing boundary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4 Calculations of ab-initio and perturbative susceptibilities

53

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

4.2

Electric Susceptibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.2.1

Ab-initio calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

4.2.2

Perturbative calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4.2.3

Results of perturbative series and ab-initio calculations for electric susceptibilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

4.3

Low-order harmonic generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

4.4

Perturbative and ab-initio calculations of low order harmonics . . . . . . . . . . . . . 67
4.4.1

4.5

Intensity dependence of low-order harmonic generation . . . . . . . . . . . . . 69

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5 Time Dependent Susceptibility

74

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

5.2

Theoretical Method

5.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

5.2.1

Time dependent susceptibility . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

5.2.2

Windowing functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

5.2.3

Segmentation of time dependent susceptibility. . . . . . . . . . . . . . . . . 83

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.3.1

Bound-bound contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

ix

5.4

5.3.2

Continuum state contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 87

5.3.3

Cross term contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

5.3.4

Populations and susceptibilities for high and low angular momentum . . . . . 90

Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

6 Conclusions

95

Bibliography98

Figures

Figure1.1

Schematic view of a filament as regions of self-focusing and defocusing . . . . . . . .

3.1

All B-Splines of order n = 3 for the knot sequence shown by full circles on the bottom

3axis. Note that multiple knots can be used in order to induce discontinuities in thebasis functions. Reproduced from [53].
3.2. . . . . . . . . . . . . . . . . . . . . . . . . 35

Energy as a function of n for different Rmax from 300 to 1000. Note the quadratic
dependence of En as a function of n that becomes more pronounced for smaller Rmax . 40

3.3

A demonstration of aliasing shown for state n = 2000, l = 0 with Rmax = 1000, with
100, 000 points on the left as an example of too few points, and 300, 000 points on
the right as an example of barely enough. For clarity, only the last 100 a.u. of the
wavefunction has been shown. Note the oscillating amplitude for the plot with too
few points. This is an example of a state that does not have enough points for an
accurate integration when finding the dipole moment. . . . . . . . . . . . . . . . . . 41

3.4

A comparison between the analytic and the calculated bound states for the groundstate (n = 1, l = 0) on the left, and n = 5,l = 1 on the right. The relative error of|Ψcalc − Ψanal |/Ψanal is included underneath each plot. These were calculated with
Rmax = 1000, 300, 000 points per grid, and Rmin = 10−6 . The plots only show a
portion of the r range of the wavefunctions so details can be seen. . . . . . . . . . . 48

xi
3.5

A comparison between the analytic and the calculated continuum states for the staten = 30,l = 1 on the left, and n = 200,l = 70 on the right. The relative error of|Ψcalc − Ψanal |/Ψanal is included underneath each plot. These were calculated with
Rmax = 1000, 300, 000 points per grid, and Rmin = 10−6 . The momentum wasfound from the calculated energy of the state, rather than the quantum number.
This compounds the error in the wavefunction, and results in a large error near theboundaries as the analytic wavefunction is not constrained to be zero due to theboundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

3.6

The relative error in the calculated energy eigenvalues compared to the analyticalenergy En = 1/(2n2 ) for l < 5. The right figure is a zoomed in version of the left,
to better show the errors at low n. The n > 16 states start to be influenced by
the bounding box, while n > 23 states have errors larger than 10%, shown by the
dot-dashed brown line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

3.7

The final wavefunctions from two different calculations using a 10 cycle sine-squared
pulse with carrier frequency 800 nm and intensity 1014 W/cm2 . Panel on the left
shows three different sizes of l absorbers for an n absorber of size 50. Good convergence is achieved, with l = 20 giving good accuracy at lower energy levels. Panel on
the right shows three different sizes of n absorbers with an l absorber of size 20. All
absorbers are identical for lower energies, however the n absorber of size 50 showssome reflection at high energies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

4.1

Susceptibility of neon as a function of intensity. . . . . . . . . . . . . . . . . . . . . . 55

4.2

The nonlinear harmonic susceptibilities of hydrogen as a function of intensity. . . . . 56

xii
4.3

a) Results of calculations for the different terms in the power series expansion forthe hydrogen atom, Equation (4.2), as a function of frequency ω of a uniform electric field. b) Same for the helium atom. Note the difference in the magnitude of
the susceptibilities for hydrogen and helium. The dashed line denotes 800 nm, thewavelength used for the majority of our calculations. . . . . . . . . . . . . . . . . . . 60

4.4a) Perturbative susceptibility coefficients χ(N ) in the power series expansion, Equation (4.2), as a function of the maximum principal quantum number nmax in thebasis set for the hydrogen atom. The radial box size was R0 = 500. The dashed line
denotes the transition from bound to continuum sets in the present basis set, while
the crosses indicate convergence of the results within 1% of the respective result at
nmax = 800. b) Same as (a), but for helium. . . . . . . . . . . . . . . . . . . . . . . . 61

4.5a) The perturbative nonlinear susceptibility terms χ(N ) I (N −1)/2 scaled by χ(3) I as afunction of the intensity of a uniform electric field at 800 nm for the hydrogen atom.
b) The same for helium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

4.6a) Comparison of results of perturbative (black line) and ab-initio calculations forthe susceptibility of hydrogen as a function of peak intensity of a 10-cycle and 30cycle laser pulse at 800 nm. b) Comparison of results of perturbative and ab-initio
calculations for the susceptibility of helium as a function of peak intensity of a 10cycle laser pulse at 800 nm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

4.7

The intensity where the relative error between the perturbative and ab-initio calculations diverge by 1% and 10%, as a function of the ground state energy of atoms with
nuclear charge 1, 1.25, and 1.5. The ab-initio calculations were done for a 10 cycle
pulse with center frequency corresponding to 800 nm. Convergence was checked in
the same manner as for previous results. . . . . . . . . . . . . . . . . . . . . . . . . . 66

xiii
4.8

Results of ab-initio numerical calculations for a low-order harmonic spectrum generated by a driver laser pulse at a central wavelength of 1600 nm, a peak intensity of
5 × 1013 W/cm2 and a pulse length of 10 cycles. The inset shows the relative error
between calculations using radial box sizes of Rmax = 500 and Rmax = 1000. . . . . . 68

4.9

(N )

Results for perturbative power series coefficients χnω for n = 1 (dashed-dotted lines),n = 3 (solid lines) and n = 5 (dotted lines) as a function of nmax . . . . . . . . . . . . 69

4.10 Results of ab-initio calculations for the integrated harmonic power (solid circles withsolid lines) for the (a) 1st, (b) 3rd and (c) 5th harmonic as a function of the peaklaser intensity of a laser pulse of 10 cycles at a wavelength of 1600 nm. The numericalresults are compared to a perturbative I n power law fit, which is matched to theab-initio results at the lowest intensity. The insets show the relative error between
ab-initio results and power law predictions with respect to the ab-initio results.

. . 704.11 Results for the ratio of higher-order terms to the lowest-order nonlinear term in theperturbative series expansion for (a) χω , (b) χ3ω and (c) χ5ω . Also shown is theratio of the sum of all higher order terms calculated with respect to the lowest-orderterm (solid lines).
5.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

The phase shift ∆Φ in radians as a function of time and transverse position in a gasof neon. On the left shown are results observed at an intensity below the ionizationthreshold, while the results on the right are at an intensity above the ionizationthreshold. Note the negative phase shift from the ionized electrons present on theright, after the probe pulse has ended (From [91]). . . . . . . . . . . . . . . . . . . . 77

5.2

The boxcar windowing function (a) in the time domain, (b) the frequency domain(log scale), and (c) a close up view of the frequency domain (log scale) with theharmonic orders denoted for a single cycle window. Note that the attenuation for
the first sidelobe is very poor (less than 10db) . . . . . . . . . . . . . . . . . . . . . . 80

xiv
5.3

Flattop window given in Equation 5.10, in the time domain (a), the frequency domain(b), and with a closeup of the frequency domain (c). Note that the wide central peakcovers 5 frequency bins. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

5.4

The Dolph-Chebyshev window with an attenuation set to 80db in the (a) time, (b)frequency and (c) a closeup of the frequency domain . . . . . . . . . . . . . . . . . . 83

5.5(a) The time dependent susceptibility of helium for Gaussian pulses with a wavelength of 800 nm and a full-width at half-maximum of 10 cycles at different intensities (100 TWcm−2 (purple), 150 TWcm−2 (red), 200 TWcm−2 (green), 250 TWcm−2(orange), and 300 TWcm−2 (blue)). The center of the driving electric field pulse isdenoted with a vertical dot-dashed line, and time is denoted in cycles of the driving field. The arrows point out the shifted peak in the susceptibility (1) and the“kink” (2) at highest intensities. Note that both the 250 TWcm−2 pulse and the
300 TWcm−2 pulse peak before the driving field. (b) A detailed plot for 300 TWcm−2 . 85

5.6a) Bound-bound contribution to the time dependent susceptibility, (hPb† µ̂Pb i (t)). It
follows the intensity envelope of the field for all intensities over the majority of thepulse. We note a small deviation at the end of the pulse for 300 TWcm−2 where theresidual population in the excited states results in a higher value of the susceptibility
at the end of the pulse as compared to the beginning of the pulse. (b) Populationof the excited states taken at the zeros of the driving field. Note the non-adiabatic
contribution to the population for driving fields stronger than 150 TWcm−2 , alongwith the residual population at the end of the pulse for 300 TWcm−2 . For both panels: 100 TWcm−2 (purple), 150 TWcm−2 (red), 200 TWcm−2 (green), 250 TWcm−2(orange), and 300 TWcm−2 (blue) pulses. . . . . . . . . . . . . . . . . . . . . . . . . 86

xv
5.7

(a) Continuum-continuum contribution to the time dependent susceptibility (hPc† µ̂Pc i (t)).
This contribution to the susceptibility is negative at high intensities (200 TWcm−2
and above). (b) The population of the continuum states taken at the zeroes of thedriving field. Note the exponential growth for 300 TWcm−2 , which starts slightly before the center of the pulse and peaks slightly after. This is followed by a decline to aresidual population at the end of the pulse. For both panels: 100 TWcm−2 (purple),
150 TWcm−2 (red), 200 TWcm−2 (green), 250 TWcm−2 (orange), and 300 TWcm−2(blue) pulses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

5.8

The population in the continuum states as a function of energy at cycle 1 for
100 TWcm−2 (purple), 150 TWcm−2 (red), 200 TWcm−2 (green), 250 TWcm−2 (orange), and 300 TWcm−2 (blue) pulses. Note the exponential decrease in population
for the 100 TWcm−2 pulse, while at the higher intensity pulses the distribution
extends to higher energies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

5.9

Cross-term contributions to the susceptibility: (100 TWcm−2 (purple), 150 TWcm−2(red), 200 TWcm−2 (green), 250 TWcm−2 (orange), and 300 TWcm−2 (blue) pulses.)
Note the large offset but subsequent small dip in the susceptibility. The dynamic
variation here is much less for high intensities than for the bound states or thecontinuum states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

5.10 Excited state population, for (a) l ≤ 1 and (b) l > 1. The lines represent resultsof calculations for 100 TWcm−2 (purple), 150 TWcm−2 (red), 200 TWcm−2 (green),
250 TWcm−2 (orange), and 300 TWcm−2 (blue) pulses.

. . . . . . . . . . . . . . . . 91

5.11 Continuum population, for (a) l ≤ 1 and (b) l > 1 on the right. The lines represent
results of calculations for 100 TWcm−2 (purple), 150 TWcm−2 (red), 200 TWcm−2(green), 250 TWcm−2 (orange), and 300 TWcm−2 (blue) pulses. Note the differencein the time dependence on the left, and on the right . . . . . . . . . . . . . . . . . . 92

xvi
5.12 The total susceptibility, split up into l ≤ 1 on the left and l > 1 on the right.
From top to bottom, the lines represent 100TWcm−2 (purple), 150TWcm−2 (red),
200TWcm−2 (green), 250TWcm−2 (orange), and 300TWcm−2 (blue) pulses. The
features we previously associated with the continuum are completely absent in thel ≤ 1 susceptibility, instead they show up in the higher l states. . . . . . . . . . . . . 94

Chapter 1

Introduction

Laser intensities have been increasing steadily since the first demonstration of the laser in1960 [1]. The first observation of second harmonic generation soon after [2] started the field of
nonlinear optics, which has exploded in the years since, and is now one of the largest areas of active
research in physics. Second harmonic generation was described using perturbation theory, denoting
the electric field of the laser as a small correction to the Coulomb electric field within the atom.
As the laser intensities available to experimentalists increased, the successes of perturbation theory
made it a widely used model for describing nonlinear optical processes.
By the late 70’s however, the discovery of ‘above-threshold ionization’ – where many photons
are absorbed, bringing electrons above the threshold of the continuum, and creating a series of peaksin the photoelectron spectrum spaced by the photon energy – had firmly marked the intensities
where perturbation theory no longer described the physics [3]. At 8 × 1012 Wcm−2 and 1 µm, theelectric fields used in these experiments were still a hundred times lower than the field strength felt
by the electron in the ground state of hydrogen, yet even this small relative field had proven toolarge for perturbation theory to account for.
Filamentation is an important nonlinear process that has been observed around this intensity,
making it an important nonlinear process to understand the transition between perturbative andnon-perturbative physics. Filamentation, the confinement and self-guiding of light propagating
through a medium due to nonlinear self-focusing, has been around almost as long as the laser
itself. The first observation of filamentation was in optically induced damage in glasses using a

2
Q-switched ruby laser by Hercher in 1964 [4], barely four years after the first demonstration of alaser by Maiman in 1960 [1]. This damage presented as a large pulverized spot, and then a “narrow
track emerging from the region of gross damage and pointing in the same direction as the incident
light” [4]. This narrow track was referred to as a filament in subsequent papers [5, 6], and thus thefield of filamentation was born.
Filamentation, as it relates to long (longer than picosecond scale - 1 ps = 10−12 s) pulse
lengths, is dominated by optical breakdown. Optical breakdown happens when the laser ionizes
enough of the atoms to cause a cascade effect: electrons are ionized and then accelerated by thelaser field and so collide with other atoms or its parent ion and ionize more electrons. The generated
plasma absorbs much of the energy of the long pulse, preventing long propagation distances. The
process of optical breakdown as the ultimate result of self-focusing and filamentation has been wellunderstood since the early 1970s [7]. However, more recently, it was observed that much shorter
pulses prevent the development of optical breakdown. In 1995, Braun et al. [8] observed that theintensities of a femtosecond (1 fs = 10−15 s) light pulses were increasing with distance: at 10 meters
downfield, the intensity was strong enough to cause damage to the mirrors in the experiment, whilea mirror placed proximal to the output of the laser would show no damage. This self-focusing
and intensity increase was the first long range demonstration of filamentation, and the first of
this behavior in air. The long range propagation was soon increased, with 50 meter propagation
distance (1996) [9], hundreds of meters (1999) [10] and kilometer distances (2004) [11] following.
Beyond the aforementioned long propagation distances, femtosecond filamentation is a hostto a slew of interesting physics. Due to nonlinear self-steepening, where the back of the pulse speeds
up in relation to the front due to the smaller index of refraction brought on by ionization, a shortlight pulse propagating through air will get shorter, as the intensity gradient at the end of the pulse
gets stronger [12]. Third harmonic generation, the shortening of the pulse from self-steepening, anda chirping of the pulse from self-phase modulation, combine to broaden the spectra of the pulse,
leading to white light generation [13]. Coherent sub-terahertz light is generated transverse to the
propagation direction in a process that is not entirely well understood [14].

3

x

z
Figure 1.1: The dynamic balance between self-focusing and ionization induced defocusing is shown
schematically. The filament is the area after the first initial focusing, while the beam is oscillating
between these two regimes. Plasma generation marks the end of the self-focusing regime and the
beginning of the defocusing regime. Figure from [15].

The physics of femtosecond filamentation (hereafter referred to as filamentation), is dominated by three processes, schematically shown in Figure 1.1: a) self-focusing, where the beamfocuses due to the nonlinear changes in the index of refraction of the medium as the beam propagates; b) ionization, where the beam reaches an intensity large enough to significantly ionize thematerial that it is traveling through; and finally c) plasma induced defocusing, where the ionizedelectrons produce a defocusing effect, causing the intensity to decrease until ionization becomes
negligible. These three effects interact in a dynamic balance, leading to the phenomenon of filamentation in air and resulting in a long lasting plasma channel behind a so-called “light bullet” [15].
Because the light pulses are so short, absorption by the resultant plasma – the dominant power
loss mode of a filament – is limited, in contrast to early picosecond and longer filaments, where theplasma absorption and subsequent optical breakdown lead to drastically shorter filaments.
This balance of self-focusing and ionization induced defocusing defines a maximum intensity
in the filament because the ionization prevents the intensity from getting too large, clamping theintensity below some material dependent value. Typical measurements have put this clamping at

4
5 × 1013 Wcm−2 for air [16]. This is above the intensity of certain strong laser field effects such as
above threshold ionization which was first measured for pulses at intensity of 8 × 1012 Wcm−2 [3],and high harmonic generation, first measured for infrared wavelengths in argon at 1013 Wcm−2 [17].
These processes are well known to be non-perturbative strong-field phenomena. Thus, we expectthat the intensity range that filamentation covers also contains the transition between perturbativeand non-perturbative effects. Consequently, the physics of filamentation is interesting because it
combines physics that we understand in a perturbative regime: harmonic generation, self-focusing,self-phase modulation, etc.; with physics that one studies in a non-perturbative setting: multiphoton ionization, high harmonic generation. Understanding better the crossover between thesetwo different regimes can give insights to improve models for filamentation. In this thesis, we arefocusing on examining this boundary.
To provide further insight into the relevance of the crossover between perturbative and nonperturbative physics we discuss the physics of self-focusing and ionization in more detail.

1.1

Self-focusing
The self-focusing of a laser is caused by an intensity dependent index of refraction for any

beam profile that monotonically decreases in intensity from the center. For example:n = n0 + f (I(r))(1.1)

where f is some monotonically increasing function of the intensity I with a zero intercept, n0 is
the linear index of refraction and r is the cylindrical radius coordinate of the beam. The makes theindex of refraction larger in the center than the outside. This larger index of refraction causes thecenter of the beam to move slower than the outside, changing the local Poynting vector towards
the area of greatest intensity, focusing the beam, and so causing a higher intensity.
Typically, this is modeled perturbatively: the function f is approximated by a power series
in I, assuming I is small (see Chapter 2 for more information):n ≈ n0 + n2 I(r) + n4 I 2 (r) + n6 I 3 (r) + . . . ,

(1.2)

5
where n2 is known as the Kerr index. At the intensities where filamentation is in the focusing regime,
the series is approximated by the first nonlinear term. This approximation is well supported byexperiments (see e.g. [18]).
In order for self-focusing to dominate, and eventually cause ionization, it needs to overpower
the defocusing inherent to any light beam – diffraction. Diffraction is fundamentally a Fourier
uncertainty relation phenomenon: the transverse momentum of the beam cannot be exactly zero,and thus a transversely localized beam will widen as it propagates in space. It can be shownthat the balance of these two phenomena, diffraction and self-focusing, will happen in any beamprofile that monotonically decreases in intensity from the center if the power is equal to the criticalpower [19]:Pcr = βλ20,
8πn0 n2(1.3)

where β is a scalar determined by the beam profile, λ0 is the central wavelength of the laser invacuum. For Gaussian beams, β is roughly 3.77 [15]. For initial powers larger than Pcr , we get selffocusing, while for initial powers less than Pcr , self-focusing will ultimately loose out to diffraction.
For argon, which has an n2 of approximately 1 × 10−19 cm2 /W at 1 atm [20], the critical power at
800 nm is roughly 1010 W.

1.2

Ionization
Self-focusing cannot lead to an intensity increase without limits. The ionization of the under-

lying medium of filamentation clamps the intensity that is possible in a filament. When the intensity
gets to a certain point – typically 1013 to 1014 Wcm−2 for air – the underlying medium gets ionizedsignificantly and the resultant plasma defocuses the beam, which decreases the intensity.
This clamping intensity is set by the process by which ionization happens. If the photon
energy is comparable to the binding energy of the underlying medium the beam is propagatingin, then ionization happens fast and the pulse quickly looses energy due to few photon absorption
processes before it can propagate very far. This corresponds to well understood perturbative ion-

6
ization. Conversely, when the binding energy is much larger than the photon energy many photons
need to be simultaneously absorbed to overcome the ionization threshold and the corresponding
probability for multiphoton absorption is very low until the intensity is large enough to ionize thematerial.
The intensity inside the filament is not easy to determine experimentally. One methodthat has been used in the past by Lange et al. [21] to determine the intensity is to introduce afilament, after propagation through air, into a cell containing argon. Inside this cell, the processof high harmonic generation takes place, and the atoms emit harmonics at odd multiple of thefundamental frequency. The largest harmonic n produced is linked to the intensity by the wellestablished relationship [22]:~nω0 = Ui + 3.2Up ,(1.4)

where Ui is the ionization energy of the underlying medium (in this case argon), and Up =e2 E 2 /4me ω02 , where ω0 = c/2πλ0 , is the ponderomotive energy of the electron in the laser field.
With this method Lange et al. found the peak intensity of a filament to be roughly 5×1013 Wcm−2 ,in agreement with numerical simulations [21].
At a high enough intensity, the atom or molecule undergoing multiphoton absorption will notabsorb just enough photons to ionize, but will also have some probability of absorbing additional
photons, resulting in the so-called continuum to continuum transitions of the involved electrons.
This behavior is popularly known as above threshold ionization [3], and its observation marked thebreakdown of perturbation theory. The breakdown intensity of the lowest non-vanishing order of
perturbation theory occurs for near optical wavelengths at 5 × 1012 Wcm−2 [23]. This is the firstindication that the ionization in filamentation has to be considered as a non-perturbative process.
High-order harmonic generation itself is another well known non-perturbative phenomena. Inhigh harmonic generation, photons are absorbed by the electron well beyond the energy needed toionize the atom, and then the electron is recombined with the parent ion, releasing a photon withenergy equivalent to the number of photons absorbed by the electron during ionization. Rather than

7
the perturbative behavior of low-order harmonic generation, where only a few photons are absorbed(typically for a total energy below the ionization threshold), and the probability of absorbing more
photons exponentially decays with the number of photons absorbed, high-order harmonic generation
generates photons with energy nω0 and (n + 1)ω0 with approximately the same probability.
The ionization of the underlying medium adds a defocusing term to the susceptibility. This
term is typically modeled as a negative correction to the index of refraction that is proportional tothe density of free electrons [15]:∆n ≈ −

ρ
,
2ρc

(1.5)

where ρ is the electron density, and ρc ≡ 0 me ω02 /e2 is the critical plasma density – the density
above which it becomes opaque. This term causes a defocusing of the beam, leading to a decrease
in intensity. This decrease in intensity eventually halts the process of ionization, which restarts thecycle, putting the beam back into the intensity regime that leads to self-focusing. The multiphoton
absorption and ionization of the pulse is the primary loss mechanism of a very short filamentation
pulse.

1.3

Higher order Kerr effect (HOKE)
Femtosecond filamentation has been around long enough that much of the basic physics has

been thought to be explored. However, in 2009, Loriot et al. challenged the traditional picture offilamentation as ionization arrested self-focusing [24]. Instead they proposed that the higher order– greater than two – Kerr terms had a negative value, leading to defocusing due to the saturationand inversion of the Kerr effect as a dominant effect, rather than ionization [25]. They did so
by extracting the higher order Kerr terms from the observation of the index of refraction using across polarization measurement that they claimed to be insensitive to the plasma contributions.
As seen in Table 1.1, from their analysis they found that the higher order terms had alternating
signs: n4 was negative, n6 was positive again, n8 negative, and n10 positive. This led them topropose that one could describe filamentation without ionization and plasma defocusing. Bejot

Table 2. Measured coefficients of the nonlinear refractive index expansion of nitrogen, oxygen, ar-gon, and air with Iinv the intensity leading to nKerrk = 0. The uncertainty corresponds to two standard
deviations of the fitted values over a set of experimental records.

Gas

N2

O2

Ar

Air

n2 (10 7 cm2 /TW)
n4 (10 8 cm4 /TW2 )n6 (10 9 cm6 /TW3 )n8 (10 11 cm8 /TW4 )n10 (10 13 cm10 /TW5 )2.2 ± 0.4
-0.16 ± 0.08
0.56 ± 0.06
-2.2 ± 0.2-

3.2 ± 0.7
-1.55 ± 0.16
1.9 ± 0.2
-10.5 ± 0.7-

2.4 ± 0.6
- 0.45 ± 0.9
0.84 ± 0.09
-4.0 ± 0.3-

Iinv (TW/cm2 )33

192.01 ± 0.19
-0.11 ± 0.31
1.6 ± 0.2
-8.6 ± 0.5
5.3 ± 0.3
26

8

26

Table 1.1: The higher order
Kerr effect values from [24] for various air components and argon.
5
0-5
Kerr||

N2
O2
Ar
then performed a simulation-20
Air of femtosecond pulse propagation that showed that one could

et al.

-10
-15
-25

-6
get defocusing, and filamentation,
without ionization using the above mentioned higher order Kerr
-30x100effect values [25].

10

20

30

40

-2

Intensity (TW.cm )Nonlinear refractive index variation of air constituents versus intensity at room temperatureThis Fig.
new3.modelof filamentation generated a significant amount of discussion. Béjot et al.
and 1 atm. (a) N2 , (b) O2 , (c) Ar, and (d) air.

and others did further theoretical studies of the consequences of higher order Kerr terms [25–28].

bound excited electrons. The high-order terms of the Kerr effect have to be taken into account
Others tried to extract the nonlinear susceptibilities from full ab-initio time dependent Schrödinger
when modelling the light propagation in air. We recall that the filamentation process in gas is
usuallydefined onby thebalancing
betweenfocusing
ionization
defocusing.
equationsolutionsa grid[29, 30],or toKerr
calculatetheandnonlinearindicesusing Thethe present
strong-fieldwork suggests that the Kerr effect contributes to the defocusing mechanism. For some gases, thefilamentation
can32].
be governedby al.
the instantaneous
effect,while
ionizationthe
plays
a minorapproximation[31,
Kolesik ettook a different
path,andthesimulated
thirdand fifthrole [15, 16]. In general, it will strongly depends on the gas and experimental conditions.
harmonic yield assuming only a third order susceptibility or the whole series of higher order Kerr
4. Conclusion
terms [33], this led them to propose a benchmark measurement for the differences between the yields
New experimental determination of the nonlinear Kerr index of refraction of N2 , O2 , and Ar
of theconstituents
third and fifth
ThisThebenchmarkwas tested,
but inconclusive
resultsare reported
have harmonics.
been performed.
time resolved
birefringence
method allows
to measure
the Kerr coefficients calibrated with the postpulse molecular alignment without the detrimental
in [34,35]. contribution.
Therefore theresultintensity,
by Loriot
et al. remained
controversial
becauseit challenged
plasmaAt highthe saturationof the electronic
Kerr effectis observed,
followed by a sign inversion above an intensity of few tens of TW/cm2 for the three gaseous
our understanding of the filamentation process, and asked questions about the validity of thecomponents. This work reports, to our knowledge, the first experimental evidence of the sign
inversionseriesof Kerr
terms at highintensity. The
result inis expected
to play a dominant role
perturbative
expansionat intensities
thatpresent
are present
filamentation.
in the self guiding of ultrashort laser pulses. In particular, the usual description of a filament
withthecalculations
plasma as main
defocusing
in consideration
the
Basic
for thevalues ofcontribution
the variousbecomes
terms ofquestionable
the perturbative
series raiseofquestions
high-order Kerr terms revealed in the present work. An illustration of their influence reported
aboutinthetheindex of ofrefraction
seriesextracted in these experiments. Using the[15]convergencedemonstratesofthe
possibility
plasma-freefilamentation.
valuesfrom the paper by Loriot et al. [24], for typical intensities of filamentation (10 TWcm−2 toAcknowledgments−2 ),wasThis work
supportedby theRégional de Bourgogne, the ANR COMOC, and the
100 TWcm
we canexaminethe Conseil
series as:th
FASTQUAST ITN Program of the 7 FP.

nAr (I = 10TWcm−2 ) ≈ 1 + 2 × 10−6 − 0.11 × 10−6 + 1.6 × 10−6 − 8.6 × 10−7 + 5.3 × 10−8 , (1.6)−2−52009; revised 15−4−3−3#110569- $15.00 USD
Apr
acceptednAr (I= 100TWcm) ≈Received
1 + 2 ×2910− 0.11 × 10 Jun+2009;
1.6 ×
10−3 15− Jun
8.62009;× 10published+ 5.321×Jul102009, (1.7)(C) 2009 OSA

3 August 2009 / Vol. 17, No. 16 / OPTICS EXPRESS 13434

9and from the increasing values we see that this perturbative series may fail to converge at intensities
relevant to filamentation. This raises the question: is the physics of filamentation accurately
described by the perturbative higher order Kerr effect?

1.4

Thesis Outline
In Chapter 2, we present a brief overview of perturbative nonlinear optics, focusing on theareas relevant to our study of the degenerate susceptibilities for comparison with our nonperturbative calculations. We start with a discussion of the principles of perturbation theory, and then fromthere extract out the response functions and the susceptibilities, being careful to ascribe physicalmeaning to the different terms where appropriate. Finally, we discuss our direct calculation of thesusceptibilities, and how these susceptibilities relate to the indexes of refraction used in the higher
order Kerr effect terms.
In Chapter 3, we introduce the numerical methods we use for the remainder of this thesis.
We discuss our method for directly solving the time dependent Schrödinger equation, using afield-free basis set for an atom in an intense laser field. These methods allow us to capture thenonpertubative physics inherent in this situation. We are specifically interested in how our method
allows for comparisons with perturbative physics by using the same basis set for both perturbativeand nonperturbative calculations.
Chapter 4 presents the comparison of the perturbative and nonperturbative calculations fora few model systems. We show that it is unlikely that the higher order perturbative series is
responsible for the saturation and defocusing seen in filamentation. We further examine the thirdand fifth harmonic perturbative series, in response to Kolesik et al. [33], examining where theab-initio calculation diverges from the perturbative calculation.
Chapter 5 discusses results from our time dependent susceptibility calculations. We areinterested in examining the susceptibility of an atom interacting with a laser field, as a functionof time within that laser field. This is applicable to filamentation studies, and high harmonic
generation phase matching. We examine the susceptibility of helium at the intensity regime where10
we previously found the divergence between the ab-initio calculations and the perturbative series
calculations, with a focus on the features that arise as we reach the threshold of this divergence.
Finally, we will conclude with a summary of the work performed, and a brief outlook towards
future work.

Chapter 2TheoryThe higher order Kerr effect (HOKE) relies on the perturbative expansion of the refractive
index in a power series, describing the self-interaction of light at a given frequency. Here, we
examine the electric susceptibility as a proxy for the index of refraction in general, and HOKE
specifically. The information in this Chapter has been discussed in many books and articles; for agood phenomenological overview, see Boyd [18], and for a strong theoretical grounding, see Butcher
and Cotter [36].

2.1

Index of refraction
The index of refraction is a macroscopic quantity relating the phase velocity v of an electricfield in a medium to that of light in a vacuum, c:
n=

v
,
c

(2.1)

For all media, the index of refraction is dependent on the frequency of the electric field of interest.
This frequency dependence arises from resonances at high frequency (typically visible light), andvibrations, rotations, ionic polarization, etc. for low frequencies (infrared, microwave and radiofrequency wavelengths).
Changes in the index of refraction that are not parallel to the direction of propagation leadto a modification of the direction of the propagating wave. This is best known as Snell’s law for an12abrupt change in the index of refraction. At the boundary of the alteration, Snell’s law gives:
nα sin(θα ) = nβ sin(θβ ),(2.2)

where nα (nβ ) is the index of refraction in the first (second) medium, and θα (θβ ) is the angle ofthe direction of propagation with respect to the normal of the interface between the two media.
The ray within the medium with the larger index will have the smaller angle. In a regime with acontinuously varying index of refraction this makes rays of light bend towards regions with a highindex of refraction, a property that is useful in gradient index lenses and is also responsible for the
self-focusing in filamentation.
For large electric field strengths, the index of refraction is no longer independent of the fieldstrength. The optical Kerr effect is due to an intensity dependent index of refraction. In theliterature, the Kerr effect is frequently limited to the first nonlinear term of the perturbation series
in intensity, however the recently discussed [24] higher order Kerr effect refers to the expansion to
all higher order perturbative terms:
n(ω) = n0 (ω) + n2 (ω)I(ω) + n4 (ω)I 2 (ω) + n6 (ω)I 3 (ω) + . . . ,(2.3)

where n is the total index of refraction, n0 is the linear index of refraction, n2 (ω) is the first orderterm of the Kerr effect (typically referred to as the Kerr effect in the literature), n4 and above are
the higher order Kerr effect, and I is the intensity.
The nonlinear index, for a positive coefficient of the dominant term, will lead to self-focusing(where Snell’s law directs the local light rays towards the axis of propagation), while for a negative
coefficient of the dominant term, de-focusing (where Snell’s law directs the local light rays awayfrom the axis of propagation) of the beam will occur. Thus, a negative n4 and n8 such as those
extracted by Loriot et al. [24] (along with a positive n2 ) and discussed in Section 1.3, would lead toa saturation of the index of refraction, as the terms with negative coefficients reach the same value
as the terms with positive coefficients at a certain field strength. As the strength of the negative
terms rises faster with intensity, this would eventually lead to an intensity regime with a negative

13
slope for the index of refraction as a function of intensity and ultimately to a defocusing of thebeam.
The index of refraction is related to the single atom susceptibility through the relation:n2 (ω) = 1 + N χ(ω),(2.4)

where χ is the single atom susceptibility of the material and is parameterized by frequency ω. N is
the number density of atoms in the material, and n is the total index of refraction for the material
at frequency ω. As shown below, the single atom susceptibility is related to the dipole moment andhence the time-dependent wavefunction of the system. It is therefore a quantity that is theoreticallyaccessible from first principles. On the other hand the index of refraction is a macroscopic quantitythat describes the interaction of light with an ensemble of atoms. Note that for the remainder ofthis thesis, we are concerned with the single atom susceptibility – where the number density of thematerial has been divided out. This is in contrast to some nonlinear optics textbooks, which areconcerned with the macroscopic susceptibility.
Transitioning between the terms of the nonlinear index of refraction and the terms of thenonlinear susceptibility is not particularly difficult. We can expand the total susceptibility andtotal index of refraction in Equation (2.4) into their respective perturbative power series:(n0 +n2 |E|2 +n4 |E|4 +n6 |E|6 +. . . )2 = 1+N χ(1) +3N χ(3) |E|2 +10N χ(5) |E|4 +35N χ(7) |E|6 +. . . ,(2.5)
where the prefactors stand for the degeneracy that we will discuss later in this Chapter. Note, forthe rest of this thesis, that the superscript in parentheticals ((n)) denotes the nth order term in
the perturbation series of a certain frequency. The electric field is related to the intensity via:
|E|2 =

I,
2n0 0 c(2.6)

where n0 is the linear part of the index of refraction, and 0 is the electrical permittivity of vacuum.
The power series can be formulated in either powers of |E|2 or powers of I, Equation (2.3) is given
in terms of I. Collecting terms with the same dependence on |E|2 gives us the first few terms of

14
the index of refraction in terms of these susceptibilities:qn0 = 1 + N χ(1) ,(2.7a)

n2 =(2.7b)

3N χ(3),4n0
10N χ(5) − 4n22
n4 =
,
8n0
35N χ(7) − 16n2 n4
.
n6 =
12n0

(2.7c)(2.7d)We are primarily interested in the role of the higher order terms, in particular how the termsof the series in Equation (2.3) behave as the intensity increases. Our analysis shows that thesusceptibility series and the index of refraction series are closely related, allowing us to use theexpansion for the susceptibility as an appropriate proxy for the behavior of the index of refraction.
We also note that any calculation of the perturbative index of refraction terms must be done throughthe susceptibilities, which makes the latter a more natural quantity for a theoretical analysis.

2.2

General properties of the electrical susceptibility
The electrical susceptibilities are used to describe the frequency dependent response – inthe form of the dipole moment – of an atom (or other form of matter) to a driving electric field.
In the simplest case, for a linear system, the susceptibility is the frequency dependent factor of
proportionality between the electric field strength and the dipole moment:p(ω) = hµ̂i = χ(ω)E(ω)(2.8)

with the driven dipole moment having an amplitude of 2 |χ(ω)| |E(ω)|, and a phase angle ofarg(χ(ω)) + arg(E(ω)).
More generally, for a susceptibility that is not confined to be linear, and an electric field E(t),with Fourier components Ej (ω), that is not confined to be a continuous wave, a total susceptibility
is related to the polarization by:
Pi (ω) = N pi (ω) = N χij (ω; E(t))Ej (ω),

(2.9)

15
where Einstein summation notation is used, and χij (ω; E(t)) is the total single atom susceptibility
tensor acting on the jth component of the electric field, giving rise to the ith component of thedipole moment and polarization. We have denoted Pi as the macroscopic polarization component
in the ith direction, and pi as the dipole moment of a single atom in the ith direction. The implicit
dependence of the susceptibility on the total time dependent form of the electric field is present
in the total non-perturbative susceptibility, but becomes an explicit dependence on the frequency
dependent components of the electric field in the perturbative susceptibility terms. For simplicity,
we consider here only with the interaction of the electric field in one, say the z, direction and thepolarization induced in the same direction, which removes the tensor indices in the above equation.
The total susceptibility χ(ω; E(t)) in Equation (2.9) can be expanded in a perturbative series
about the small quantity E(ω) as:χ(ω; E(t))E(ω1 ) = χ(1) (ω; ω1 )E(ω1 ) + χ(2) (ω; ω1 , ω2 )E(ω1 )E(ω2 )+(2.10)

χ(3) (ω; ω1 , ω2 , ω3 )E(ω1 )E(ω2 )E(ω3 ) + . . . ,where each susceptibility term χ(n) is parameterized in terms of the frequency ω of the dipoleP
moment it drives and the electric field frequencies ωj that contribute to that term, where ω = j ωj .
Using Equation (2.10), and Equation (2.9), we can write down the perturbative terms for the

expectation value of the dipole moment for some frequency ω as a sum of the contributions due toelectric field components at different frequencies:hµ̂(n) (ω)i =X

χ(n) (ω =ω1 ...ωnω=ω1 +···+ωn

n
X
i=1

ωi ; ω1 , . . . , ωn )E(ω1 ) · · · E(ωn ),(2.11)

where µ̂(n) is the nth order term in the perturbative expansion of the dipole moment. This gives
us a way to obtain the susceptibility terms, given the Fourier decomposition of the perturbative
dipole moment terms, by collecting terms.
The perturbative expansion of the dipole moment, via the expansion of the susceptibilities,
is not always desired. Alternatively, the total susceptibility can be found from ab-initio calculations for a given dipole moment and a given electric field as a function of time. By rearranging

16
Equation (2.9) we getχ(ω; E(t)) =p(ω)
F [hp(t)i] (ω)=,E(ω)
F [E(t)] (ω)(2.12)

where we make explicit that this susceptibility depends on the time-dependent electric field, andF denotes the Fourier transform of its argument, allowing us to convert from an arbitrary timedependent dipole moment to the frequency dependent susceptibility. This is how we extract out
the susceptibility from our ab-initio calculations.

2.3

Perturbation Theory

2.3.1

Brief review of perturbation theory

We derive the perturbative susceptibilities via the Fourier transform of the time dependent
terms of the perturbative dipole moment, i.e. the susceptibility is extracted from Equation (2.11).
The perturbative dipole moment terms are found using the time-dependent interaction with theelectric field (HI (t)) as a small perturbation to the field-free Hamiltonian H0 . For the expansion
we write the corresponding Schrödinger equation as:(H0 + λHI (t))Ψ = −i∂t Ψ,(2.13)

where Ψ represents the wavefunction, λ is a parameter that will be used to denote the orders in
the perturbation series, and will be set to 1 at the end of the calculation. If Ψ is expanded into a
power series in λ:Ψ = ψ (0) + λψ (1) + λ2 ψ (2) + . . . ,(2.14)

where ψ (n) is the nth order term in the expansion of the wavefunction Ψ, an application of Equation(2.14) along with Equation (2.13) allows for the collection of terms in powers of λ and the derivationof the recurrence relation (after setting λ to 1) :H0 ψ (0) = −i∂t ψ (0)(2.15)

H0 ψ (n) + HI ψ (n−1) = −i∂t ψ (n) .

(2.16)

17
By decomposing the perturbative terms of the wavefunction in terms of the eigenstates φmof H0 , where ωm = ~Em is the frequency associated with that eigenstate:ψ (n) =

X

−iωm tc(n)m (t)φm e(2.17)

mwe obtain a series of recurrent differential equations:ċ(0)m (t) = 0−1ċ(n)m (t) = (i~)(2.18a)X (n−1)ck(t) hφm | HI (t) |φk i eiωmk t

for n ≥ 0,k

(2.18b)which can be recast into integral form as:−1c(n)m (t) = (i~)

Z t

−∞

dτ

X
k

(n−1)

hφm | HI (τ ) |φk i ck

(τ )eiωmk τ(2.19)(0)with cm is the time independent coefficient that defines the initial conditions of the problem.

2.3.2

Dipole moment and single atom response functions

Up until now, the form of the time-dependent perturbation HI (t) has not been made explicit.
We are concerned with the interaction of an atom with an oscillating electric field in this thesis,
so we proceed to write the time-dependent dipole moment of the system as the response of the
field-free system to an electric field.
For simplicity, we make three assumptions here, which we will discuss in more depth in Section 3.1, as they are of significant importance when solving the ab-initio time-dependent Schrödinger
equation. First, we choose linear polarization along the ẑ direction for the electric field interaction. This reduces the complexity of the following equations, without losing any relevant meaning
towards our discussion, and matches the scenarios considered in the ab-initio calculations. Second, we assume a single active electron potential, which are used in the ab-initio calculations, forobtaining the basis functions (see Section 3.1.1 for more information on the single active electron
models). Finally, we use the length gauge for the interaction with the electric field. This simplifies
the relationship between the dipole moment and the interaction Hamiltonian, and again matches

18
the ab-initio calculations performed in this thesis. The interaction Hamiltonian is thus given by:
HI = −µ · E(t) = −µE(t)(2.20)and the perturbative expansion of the time dependent dipole moment yields:
hµ̂i = hΨ(t)| µ̂ |Ψ(t)i(2.21)hµ̂(0) i}|{z= hψ (0) | µ̂ |ψ (0) i +hµ̂(1) iz}|{hψ (1) | µ̂ |ψ (0) i + hψ (0) | µ̂ |ψ (1) i +hµ̂(2) iz}|{hψ (2) | µ̂ |ψ (0) i + hψ (1) | µ̂ |ψ (1) i + hψ (0) | µ̂ |ψ (2) i + . . .
The zeroth order term hµ̂(0) i is the permanent electric dipole moment, which is zero for all centrosymmetric systems including the atomic systems we are studying. The first order correctionhµ̂(1) i is the non-zero linear response of the system to the applied electric field. Application of
the recurrent relations in Equation (2.19), and the assumption that the system we are studying is
(0)

(0)

(0)

prepared in the eigenstate |gi, i.e. ψ (0) = cg φg e−iωg t with cg = 1 and cm6=g = 0, gives the linear
term of the dipole moment:

hµ

(1)

(t)i =

Z t

−∞

z X

1
dτ i~

m6=g

R(1)(t−τ )}|{−iωmg (t−τ )µ̂gm µ̂mg e
E(τ ) + c.c.,(2.22)

where µ̂ij = hφi | µ̂ |φj i is the matrix element of the dipole moment operator corresponding to thestates denoted by the labels i and j, and ωij = ωi − ωj is the energy difference between thesestates. Equation (2.22) is the time-dependent convolution of the electric-field response of the atomwith the electric field itself. The integration limits ensure that R(1) (t − τ ) is never evaluated fornegative arguments and hence causality is guaranteed. Further, the response is constrained to bea function of the difference in time between the current time t and the time τ of the interaction
with the electric field, so that the polarization is invariant with respect to a global time offset. The
response is a weighted sum of the phase factors of the various eigenstates of the system.

19
Using
∞
X

E(t) =

E(ωj )e−iωj t ,(2.23)

ωj =−∞

we can evaluate the integrals by substituting Equation (2.23) in Equation (2.22):hµ(1) (t)i =zXXχ(1) (ω1 )ω1 m6=g}|{|µ̂gm |2
E0 e−iω1 t .
ωmg − ω1

(2.24)The term is proportional to E(t), as is expected for linear polarization. Note, there is a simple
pole at ω0 = ωmg , related to resonance behavior, when the driving field has the same frequency as
the atomic transition. However, Equation (2.24) neglects the lifetimes of the atomic states. Whentaken into account, the lifetime of a state will broaden the transition and remove the singularity.
This effect is frequently approximated by replacing a real ωmg by ω̃mg = ωmg + iΓmg , where Γmg
represents the spontaneous decay of the population of the state back to the ground state. This
shifts the pole into the complex plane, and therefore removes the singularity. Here, we are only
concerned with transitions far from resonance so we do not have to consider this case. Further
details on this transformation can be found in [18, 23, 36].
Examining the second order contribution to the polarization, via application of the recurrence
relation in Equation (2.19), one can start to see some patterns and the corresponding physicalmeaning of these terms. Keeping the initial condition the same one getshihµ(2) (t)i = hψ (2) (t)| µ̂ |ψ (0) (t)i + c.c. + hψ (1) (t)| µ̂ |ψ (1) (t)i(2.25a)


Z t
Z τ
X
−2
0 −iωlg (τ 0 −t) −iωml (τ −t) ∗∗ 0= (i~)dτdτ ee
E (τ )E (τ ) + c.c.
µ̂gl µ̂lm µ̂mg
−∞

m6=g
l6=g

−

Z t

−∞

−∞

∗

E (τ )e−iωlg (τ −t)

dτ

Z t

−∞

iωmg (τ 0 −t)e0E(τ ) dτ0.
(2.25b)Note that the response is symmetric with respect to an interchange of the time variables. In particular, in the case of an electric field of linear polarization, as considered here, such an interchange
cannot change the physics of the situation. This property is known as the intrisic permutation
symmetry.

20
The hµ(2) i term consists of an odd number of dipole transition matrix elements. Therefore,
the term is zero for centrosymmetric systems because the dipole selection rules for atoms require l tochange by ±1. Thus, an even number of dipole transitions is necessary to result in a non-zero term.
Nonetheless, an examination of the form of this dipole term provides us with an understanding forthe higher order responses.
The first integral in Equation (2.25b) can be related to two photon sum and difference processes such as second harmonic generation and sum frequency generation. Physically, the integration over τ 0 establishes a phase difference between the populations in the l state and the g state thatevolves with time. The second interaction with the field, which is represented by the integration
over τ , modifies this phase difference as the field again interacts with the states. This transition isenhanced if there is an intermediate state m that is close to a frequency of the laser field. Note that
this is a property of the second order perturbative wavefunction, which means that this behavior– the sum of two photons – shows up as well in third order and higher order dipole moment terms
as an intermediate process. The second integral represents two different states being excited insingle-photon processes.
We can evaluate the integrals by substituting Equation (2.23) in Equation (2.25b):hµ(2)(t)i =X X

ω1 ω2 lm−+µ̂gl µ̂lm µ̂mg
E(ω1 )e−iω1 t E(ω2 )e−iω2 t(ωlg − ω1 − ω2 )(ωmg − ω1 )µ̂gl µ̂lm µ̂mg
E(ω1 )e−iω1 t E ∗ (ω2 )eiω2 t(ωlg + ω2 )(ωmg − ω1 )

µ̂gl µ̂lm µ̂mg
E ∗ (ω1 )eiω1 t E ∗ (ω2 )eiω2 t .
(ωlg + ω2 )(ωmg + ω2 + ω1 )(2.26)The two photon sum process is shown in the first and last terms. The complex exponentials describe
a wave with frequency ω1 + ω2 , while the coefficient is singular for states with frequency ω1 + ω2 ,denoting a maximal transition amplitude. This behavior is enhanced if either ω1 or ω2 is resonantwith an intermediate state. As mentioned before, due to the lack of modeling of the lifetime of thestate, electric field components at on-resonant frequencies will produce unphysical dipole moment
terms. The second term describes the difference frequency contribution, where the polarization is

21
at the difference frequency of the two electric field terms. Together, for a two color electric field,
the dipole moment will have terms proportional to four different frequency components:− ω1 + ω2

ω1 − ω2

ω1 + ω2− ω1 − ω2 ,(2.27)

leading to real polarizations with two frequencies: ω1 + ω2 and ω1 − ω2 . When ω1 = ω2 a DC
contribution and a second harmonic contribution are produced.
Real electric fields have real Fourier coefficients E(ωj ) such that E(ωj ) = E ∗ (−ωj ). Thispermits positive and negative frequencies to be swapped without changing the total sum. This
freedom allows the electric field and exponential time dependence of Equation (2.26) to be extracted
from the sum over the dipole terms:"
XX
hµ(2) (t)i =
PI
ω1 ω2

lm

+
+

µ̂gl µ̂lm µ̂mg
(ωlg − ω1 − ω2 )(ωmg − ω1 )µ̂gl µ̂lm µ̂mg
(ωlg + ω2 )(ωmg − ω1 )

µ̂gl µ̂lm µ̂mg
E(ω1 )E(ω2 )e−i(ω1 +ω2 )t ,(ωlg + ω2 )(ωmg + ω2 + ω1 )(2.28)

where the part in square brackets is the susceptibility χ(2) (ω3 = ω1 + ω2 ), which can be seen in
Equation (2.11). Note that due to intrinsic permutation symmetry, the order of the electric fieldterms E(ω1 )E(ω2 ) does not matter. The symmetry requires that we average over interchange ofω1 and ω2 , which is represented by the operator PI . This susceptibility tensor is thus an 2nd ordersymmetric tensor with a term for every ω1 and ω2 .
A comparison between the first and second order terms of the dipole polarization shows
how the pattern of the third and higher order perturbative terms evolves. The nth term in theperturbative series will have n + 1 terms in the sum:hµ(n)(t)i =nX

m=0hψ (m) | µ̂ |ψ (n−m) i .

(2.29)The perturbative wavefunction of nth order will contribute an n photon process to the series, alongwith all intermediate state processes (processes involving 1 . . . n − 1 photons):
|ψ

(n)

−n

i = (i~)X

jn 6=j0n
X Yµjm jm−1 E(ωm )e−iωm t
P···|ψj0 i ,ωjm j0 − ni=m ωij1 6=j0 m=1(2.30)

22
where jm denotes a field-free eigenstate with j0 denoting the ground (or initial) state.
The nth order perturbative susceptibility term can be found via the same method used forχ(2) above. This includes an interchange of ωi with −ωi where appropriate to consolidate theelectric field terms, and then an extraction of the specific susceptibilities desired.

2.4

Summary
Here we have given an overview of perturbation theory as it applies to the dipole moment

and the susceptibility. Since we are interested in adding to the ongoing discussion of the higher
order Kerr effect, which was presented as a function of the index of refraction, we discussed the use
and relevance of the perturbative susceptibility as a more theoretically accessible quantity.
For the discussion of the relevance of the perturbative approach we derived the first few
terms of the dipole moment under a set of certain assumptions – such as a single active electron
potential, linear polarization, etc – that we will use throughout this thesis. Each term was discussedwith respect to its physical significance. Finally, we discussed the method that we will use in
our calculation and analysis of the susceptibilities. Results of this analysis will be presented in
Chapter 4.

Chapter 3

Numerical Methods

Solving the time dependent Schrödinger equation (TDSE) for atoms in strong electric fields
is challenging, since strong-field ionization and other effects are complex highly nonlinear processes.
An electron interacting with a powerful enough laser can simultaneously cause a number of related
phenomena including: over the barrier ionization, bound-bound excitation, high harmonic generation, above threshold ionization, dispersion of the ionized electron wavepacket and others. Thesephenomena interact in unexpected ways, leading to incredibly rich and complex physics. This ispart of the excitement of working in this field.
These complex interactions are also what makes it difficult to perform calculations. The
wavefunction of an electron interacting with a strong electric field is spread out over a muchgreater volume than a typical atom in the ground state. This volume must be represented in any
calculation that aims to accurately model this interaction.
Furthermore, the Schrödinger equation of an atom in a time dependent strong electric field
does not have an exact analytical solution, and the large number of coupled dimensions – three foreach electron – combined with the large interaction volume, requires far too much information for
a complete computational solution using current computers. Due to the difficulty of the problem,
any realistic solution is a balance between three competing factors:• Numerical and physical accuracy of the method. When dealing with differential equations
without known analytical solutions, we can use computers to find a numerical approximation to the exact solution to that differential equation. How closely these approximations

24
are to the exact solution is the numerical accuracy of the method. More accurate methods
can be more difficult to calculate, either requiring much more memory or computational
power. How accurate the numerical results are must be balanced against the accuracy ofthe physical model with which one is trying to solve the problem. The accuracy requiredby the physics that we want to study determines the physical model that we use and thenumerical accuracy we require.
• Information extracted from the method, such as physical observables, or quantum mechanical properties we wish to examine. Certain observables are easier to calculate in somemethods than in others.
• Time to calculate.
These limitations have become significantly less severe as computational power and algorithmic advances have developed over the last few decades. These advances have allowed theorists tocontinually push the boundaries of accuracy, and allowed us to examine larger systems or more
subtle observables.
In our case, we examine the perturbative and the ab-initio solutions to the time dependent
Schrödinger equation for an atom in a time dependent electric field. In order to study the limitsof the perturbative approach by comparing the respective results to those of ab-initio solutions,
we aim to perform both calculations within the same theoretical framework. To this end, we
have chosen the representation of the wavefunction in terms of the field-free basis states. Since
we are specifically interested in establishing the limits of the perturbative HOKE expansion, thisdetermines the susceptibility as the observable we need to evaluate. The susceptibility is determined
via the dipole moment (see section 2.2), which can be well calculated using the basis states. In this
Chapter we will present the numerical basis state method used to perform the perturbative andab-initio calculations of the atomic susceptibilities for the comparison.

25

3.1

Overview of the method
The Schrödinger equation for a single atom in an oscillating electric field within the dipole

approximation and using the length gauge can be written as (in atomic units):


 X
N  2
X
pi
1
1

 Ψ(r1 , . . . , rN , t) = −i∂t Ψ(r1 , . . . , rN , t),− + ri · E(t) +2ri|ri − rj |
i=1

(3.1)

i6=j

where pi = −i∇i is the usual momentum operator for electron i, with position coordinate ri ,interacting with a time dependent electric field vector E(t) and Ψ(r1 , . . . , rN , t) is the wavefunction.
The solution of the general Schrödinger equation for atoms with more than two electrons is
too computationally intensive for strong-field calculations, and must be reduced before it can besolved in a reasonable amount of time. We make a number of approximations to the multi-electron
Schrödinger equation in order to restrict the calculation time and to make efficient the calculations
that we are interested in performing.
The approximation with the largest effect on computation time is the single active electron
approximation. For hydrogenic systems, this is an excellent approximation: the outer electron ismore weakly bound and will thus react to the external electric field more freely than the innerelectrons. A good example is Lithium, with two electrons in the inner s-orbital, and a single
electron in the outer s-orbital. The outer s-orbital electron exhibits a much weaker (approximately
11 eV less tightly bound) binding energy than the inner electrons. However, even for atoms withmore than one electron in the valence shell, such as the noble gas atoms, the single active electron
model has been proven reasonably accurate in the past in the case of studies in strong fields. For
example, published comparisons between theoretical and experimental studies of the multiphoton
ionization of xeon and ATI spectra of argon have shown good agreement using single active electron
models [37, 38]. This will be covered in more depth in subsection 3.1.1.
Next, we consider only spherically symmetric potentials in order to work within the separable
spherical coordinate system. Spherically symmetric potentials include those atoms with a complete
outer electron shell such as noble gases – which are popular for experimental studies of strong
field processes – and atoms with an outer s orbital shell, such as alkali atoms and helium. Using26
spherically symmetric potentials allows us to separate the angular part of the field-free Hamiltonian– with the well known spherical harmonics as solutions – from the radial part. This requirement does
limit us from being able to model atoms with incomplete valence p-orbitals or those with higher
orbital angular momentum due to the symmetry breaking the corresponding potentials exhibit
under the influence of an external electric field. Such potentials have coupled terms between theirradial and angular parts, which significantly increases the complexity of the field-free Hamiltonian.
We select a linear polarization direction for the interacting electric field (denoting it ẑ). This
reduces the quantum numbers (and thus the dimensionality) that change during the interaction of
the atom with the field via the dipole selection rules to the principal quantum number n and theorbital angular momentum quantum number l. The angular momentum quantum number in the z
direction (denoted by m) is fixed by the choice of the initial state (typically we set m = 0).
Finally, since we will represent the radial part of the wavefunction on a grid, we need torestrict its size. To this end, we use an infinite spherical well with radius R max , beyond which thepotential is infinite. In practice, we denote this by adding a BC(r) term in the Hamiltonian which
fixes the value of the wavefunction at Rmax to zero, enforcing the boundary condition.
These approximations give us the following form of the TDSE:
 2

p
e
+ V (r) + zE(t) Ψ(r, θ, t) = i∂t Ψ(r, θ, t),
2

(3.2)

where now Ve (r) = V (r) + BC(r) is the spherically symmetric single active electron potential

plus boundary conditions, zE(t) is the electric field interaction term that is now limited to alinear polarization with z = r cos θ, and Ψ(r, θ, t) is the wavefunction parameterized in the three
independent coordinates.

3.1.1

Single active electron models

The single active electron (SAE) approximation consists of two parts: the wavefunction of thesingle propagated (active) electron, and the potential in which it is propagated. Because the SAE
approximation assumes that all electrons other than the active electron are frozen in their groundTarget
He

Zc
1.0

a1
1.231

c1
0.662

a2
-1.325

c2
1.236

a3
-0.231

c3
0.480

27

Table 3.1: The parameters from [43] for the helium atom.

state configuration, the potential in which the active electron is propagated includes contributions
from both the bare nuclear charge as well as the frozen electrons.
A SAE potential can be determined by finding the full multi-electron ground state for theatom using a method such as Hartree-Fock [39], or density functional theory [40, 41], and thenremoving the active electron. The resulting multi-electron wavefunction describes the inactive
electrons. The complete potential in which the active electron is propagated results from theelectrostatic repulsion of the inactive electrons together with the electrostatic attraction of the core
nucleus.
For spherically symmetric potentials, potential energy surfaces found in this way are typically [42, 43] fit to a series of radial Yukawa potentials:
V (r) = −Zbare X an e−cn r
+
,
r
r
n

(3.3)

where Zbare is the asymptotic bare charge that the electron will see, and an and cn are fits to thecalculated potential. In this thesis we have performed calculations for the helium atom using thesingle active electron potential from Tong and Lin [43], shown in Table 3.1.

3.1.2

Separation of variables

We want to consider a model system that can be used for ab-initio time propagation as wellas time independent perturbative calculations. This requires the separation of the time dependent
Hamiltonian into a time independent Hamiltonian plus the time dependent field interaction, zE(t).
The time independent part of the Hamiltonian is then diagonalized, giving us a time independentbasis set that we can use for the perturbative calculations, and that we can use as the elements inwhich we propagate the total time dependent Hamiltonian.

28
We find the basis state functions of the Schrödinger equation for the time independent Hamiltonian given by:




−∇2
e
H0 Ψ(r, t) =+ V (r) Ψ(r, t) = −∂t Ψ(r, t),
2

(3.4)

where Ve (r) is the spherically symmetric potential plus boundary conditions. In the sphericalcoordinate system, this Schrödinger equation is separable, allowing us to write down the solution
as a product of functions of different coordinates [44]:Φ(r, t) = R(r)Y (r̂)Ξ(t),(3.5)

where R(r) is a function dependent only on the radial distance between the electron coordinate
and the nuclear coordinate, Y (r̂) depends on the angular coordinate, and Ξ(t) specifies on the timecomponent of the wavefunction. This separation of variables is only possible because our total
potential Ve (r) depends only on the radial coordinate.

From this, we get three separate equations relating the three components:∂t Ξ(t)
Ξ(t)
h
i
1 2 2
l(l + 1) =r ∇r R(r) − 2r2 Ve (r) − E
R(r)
1 2 2
l(l + 1) = −r ∇r̂ Y (r̂)
Y (r̂)
E = −i

(3.6)
(3.7)
(3.8)

where E is the energy eigenvalue, l is the total angular momentum quantum number, ∇2r is theradial part of the Laplacian and ∇2r̂ is the angular part. This reduces the field-free problem tothree eigenvalue problems, two of which have well known solutions. The solution to Equation (3.6)
is the exponential time dependence of a state:
Ξ(t) = e−iEt ,(3.9)

while the solution to Equation (3.8) is given by the spherical harmonics:
Yl,m (r̂) = N eimφ Pl,m (cos θ),(3.10)

where N is a normalization factor, l and m are the total angular momentum quantum number andprojected angular momentum quantum number respectively, and Pl,m are the associated Legendre

29
polynomials. The radial equation (Equation (3.7)) is now defined in terms of the sphericallysymmetric potential and boundary conditions.

3.2

Numerical basis state method
The time independent calculation of the eigenstates is simplified by using the spherical har-monics (as discussed in Section 3.1.2). This reduces the calculation to merely finding the radial
part of the eigensystem from Equation (3.7):∂rrl(l + 1)e−
+ V (r) +Rn,l (r) = En,l Rn,l (r).
2
2r2

(3.11)The rest of this section is related to the solution of this eigenvalue equation. We start with a general
discussion of representations of differential equations in section 3.2.1. In section 3.2.2 we give a
brief overview of our numerical method for finding the eigenstates. This leads us into section 3.2.3,
where we discuss the free parameters that are useful for the basis. Section 3.2.4 provides details of
the algorithm that we use to solve the eigensystem, along with a comparison with analytic results(section 3.2.4.4).

3.2.1

Representations of differential equations

Much work has gone into finding different strategies for solving the time-independent Schrödinger
equation [45–47]. Like all problems in applied math, different methods have different strengths andweaknesses. Here we will discuss some of the more common techniques for solving the time independent Schrödinger equation. The general form of the Schrödinger equation for our model iswritten as:p2e+ V (r) Ψ(r) = EΨ(r),
2m
|
{z
}



(3.12)

H0

The first decision that must be made when solving any differential equation is the numericalrepresentation of the solution, and consequently, the differential operator. The wavefunction of anelectron is represented as a complex valued function defined over the three dimensional real-valued

30
spatial vector. However, when trying to solve it numerically on machines with finite memory andcomputational power, we must make choices in order to describe such a function. These choices
determine the form of our solution. We are primarily interested in this section with the representation of the eigenfunctions of the time independent Hamiltonian, and will focus our discussion on
methods to solve such a system. We will discuss the propagation of the wavefunction in time insection 3.3.
There are many choices one can make on how to represent the wavefunction and the Hamiltonian. These methods generally fall into two different categories: approximations via the discretization of space and approximations via the truncation of series. Below we will give a brief overview
of these general methods. In order to simplify our analysis, we restrict ourselves to the eigensystem
consisting of the 1-D (radial) time independent Hamiltonian, and the corresponding 1-D radial partof the wavefunction. This matches the problem we are looking to solve.

3.2.1.1

Discretization of space

One way to solve a time-independent differential equation is to discretize the variables [48].
In the case of the Schrödinger equation, this means discretizing the space in which an electron (ora center of mass) is represented. For example, for a Cartesian coordinate system, one defines agrid of points on which the solution is defined, while between these points the solution is inferred.
Typically, an evenly spaced grid is used for simplicity, however this is by no means a requirement.
For an evenly spaced grid, the distance between adjacent points – the grid spacing – we denote by∆x, and the nth point on the grid is given as xn . The total number of points is represented by
nmax .
Derivatives on such a discrete grid are called finite differences, and thus the methods thatsolve differential equations on a finite grid are called finite difference methods. Finite difference
methods are perhaps the simplest of the various methods to solve the Schrödinger equation. In
position space, the momentum operator is represented by a derivative p = −i∇, and thus the
kinetic energy term is proportional to the second derivative. The first derivative of a function f (x)31
between two adjacent points of a grid is then approximated by:∂x f (x) =f (x + ∆x) − f (x)+ O(∆x2 ).
∆x(3.13)

Equation (3.13) is referred to as the forward difference, because it approximates the derivative at
x using the value of the function at x + ∆x and x. The backward difference is easily inferred fromthe forward difference.
Combining the forward difference and backward difference, one gets the central difference fora second order derivative:∂xx f (x) =f (x − ∆x) − 2f (x) + f (x + ∆x)+ O(∆x2 ).
∆x2(3.14)

Using this, the time independent Schrödinger equation in one dimension can be written as:−HΨ(x) = EΨ(x)(3.15)1 Ψ(xn−1 ) − 2Ψ(xn ) + Ψ(xn+1 )+ V (xn )Ψ(xn ) = EΨ(xn ).
2∆x2(3.16)The boundary conditions are set by fixing the values of Ψ(x0 ) and Ψ(xmax ).
The Hamiltonian H0 in Equation (3.16) above can be represented as a matrix:
1
 m∆x2 + V (x1 )
H0 = 







1
− 2m∆x
2

0

···

0

1
− 2m∆x2

1+ V (x2 )m∆x2

1− 2m∆x
2

0

···

0

0
..
.

1
− 2m∆x2

1+ V (x3 )m∆x2

1− 2m∆x
2

0

···

.

.

..

..

..

.

..

.

..

.





,





(3.17)

where the boundary conditions have been fixed via neglecting to add the respective points to thematrix, which is equivalent to setting them to zero. Equation (3.16) can be rewritten as a matrix
equation:
Ho Ψ = EΨ

(3.18)

where the wavefunction is now represented as an n dimensional vector. This formulation allows
one to use matrix eigenfunction solvers.

32
The finite difference method in its simplest form is very easy to implement and reason about.
For low dimensional calculations, it is compact in its representation, and with a small enough grid,
it usually has acceptable errors. Unfortunately, the finite difference method is difficult to adapt
to non-uniform grids [48, 49], and with a large number of dimensions, its requirement for a small
grid-spacing quickly becomes computationally unfeasible. Finite difference methods are frequently
used as well to solve the time-dependent Schrödinger equation, using a finite difference to discretize
the time component as well as the space components.

3.2.1.2

Finite series representations

Another method to represent the wavefunction and Hamiltonian is as an abbreviated series.
Basic calculus allows one to approximate any well behaved arbitrary function as some finite numberof terms in a power series expansion such that the desired accuracy is achieved. More generally, one
can represent any well behaved arbitrary function f (x), in this case the wavefunction, as a truncatedseries in any complete basis S with basis elements sn (x) (for simplicity, we have restricted ourselves
to a discussion in one dimension):f (x) =

X

cn sn (x),(3.19)nwhere cn is a coefficient that represents the weight of the individual functions sn (x) within the sum.
This allows one to potentially represent f (x) very accurately in a small amount of memory, given
one chooses an appropriate basis set S, such that the number of basis states needed for convergenceis small.
This scheme can be expanded to time dependent problems by defining a time dependent
coefficient to each basis element:f (x, t) =X

cn (t)sn (x).

(3.20)n

A linear operator Ô can be represented within S as a matrix with weighted transitions
between basis elements sn and sn0 in some domain R:
Ôn0 ,n =ZRs†n0 (x)O(x)sn (x) dx.

(3.21)

33
Consequently, the Hamiltonian within the Schrödinger equation can be represented as:Z0s†n0 (x)H(x)sn (x) dx,
Ĥn ,n =(3.22)

R

for a basis set S which consists of orthonormal elements sn . When S is not orthonormal, then s†nis no longer the complex conjugate of sn , and must be constructed out of the overlap matrix (formore information, see [50, 51]).
From a computational perspective, a low memory footprint of the operator is important.
There are two properties of an operator that can lead to a low memory footprint: the total numberof states needed to accurately represent the wavefunction and the sparseness of the operators acting
on such a wavefunction. If one knows, a priori, that certain elements will be zero within the matrix,
then such elements do not need to be represented, and storage can be reduced.
There is an upfront cost to using a finite series representation. The Hamiltonian must be
found within this representation using Equation (3.22) before one can do any further calculations– such as propagation – in this basis. Typically, the Hamiltonian in this form can be saved andused for multiple further calculations. When using the field-free basis, the Hamiltonian must berecalculated for every combination of atomic species and grid that one wishes to use.
A few basis sets are particularly common for ab-initio solutions of the time dependent
Schrödinger equation in problems in atomic, molecular and optical physics. Each basis set has
benefits and drawbacks. In general, the sparsity of the Hamiltonian within the representation –the number of non-zero matrix elements compared to the size of the matrix – and the number of
elements of the basis set S necessary to approximate the wavefunction determine the memory andcomputation required for a calculation.
A common expansion used is the Fourier series [52], where an arbitrary complex function isdecomposed as a sum of complex exponentials over the interval of convergence [x0 , x0 + P ]:f (x) =∞
X

2πncn ei P x .

(3.23)n=−∞

By truncating the series to a finite number of exponentials, the sum represents an approximationof the function. For an atomic, field-free potential, the Fourier representation of the Hamiltonian

34
is very dense: the matrix element between any two frequencies is nonzero, for a typical atomic
potential proportional to −1/x. Furthermore, the representation of the ground state of a hydrogen
atom requires a large number of frequency components to reach convergence to sufficient accuracy.
On the other hand, for a free electron without any potential (H = p2 /2), the Hamiltonian is
diagonal in the Fourier representation. The Fourier basis functions are eigenfunctions of the kinetic
energy term of the Hamiltonian, allowing for an efficient description of the free electron problem
within this basis set.
We have already used the spherical harmonics representation in section 3.1.2. For sphericallysymmetric problems, the angular parts of the wavefunction are easily represented in terms ofspherical harmonics. For atoms interacting with an electric field, the dipole transition rules (l0 = l±1
for the total angular momentum quantum number l and m − 1 < m0 < m + 1 for the magnetic
quantum number m if necessary) give a very sparse matrix representation of the dipole transition
operator for single active electron potentials. Since the spherical harmonics are eigenstates of theangular part of the single active electron field-free Hamiltonian, the matrix representing this partof the field-free Hamiltonian is diagonal.
Another widely used method to represent the radial part of a wavefunction are B-Splines [53,54]. B-Splines are a family of piece wise polynomials defined via the recurrence relation:
1 if ti ≤ x < ti+1
Bi,1 (x) =0otherwise
Bi,n (x) =x − titi+n − xBi,n−1 (x) +Bi+1,n−1 (x),ti+n−1 − titi+n − ti+1

(3.24)(3.25)

where Bi,n (x) is the B-Spline of order n defined between the “knot”s ti and ti+k . These knots canbe at arbitrary points within the space that is being used, and multiple knots at the same point canbe used to fulfill boundary conditions and allow discontinuities to be represented. This behavior
can be seen in Figure 3.1, where each set of lines represents a basis function of order 3 that iscentered at different knot points. The locations of the knots are shown at the bottom axis. Noticethat at multiple knots (at the far right and far left) a B-Spline may have discontinuous derivatives.

35This is useful for representing Neumann boundary conditions and discontinuous potentials.
For a B-Spline basis the accuracy improves as the density of the knots increases, similar to theincreased accuracy when decreasing grid point spacing in finite difference schemes. B-Splines arefunctions with limited support – they are only non-zero within a small number of knots – leading
to a sparse representation of the radial part of the Hamiltonian. B-Splines of nth rank are nthorder polynomials, which have a well defined nth derivative when defined over n knots with uniquelocations. Thus, the transition elements of the Hamiltonian can be obtained analytically for ananalytical potential. Unfortunately, B-Splines do not form an orthonormal basis, which means thatan overlap matrix must be calculated in order to solve the time-independent Schrodinger equation,1820
as mentioned
above.

H Bachau et al

x
0.0
1

1.0

2.0

3.0

4.0

5.0

B3(x)
0.8

B(x)

0.6

0.4

0.2

0

Figure of1.order
Then fullof ordershownk = by3 fullrelativeto theknot bottom
sequence
Figure 3.1: All B-Splines= 3 setfor oftheB-splines
knot sequencecircleson the{0, 0, 0, 1, 2, 3, 4, 5, 5, 5}. Knots are represented by full circles.
axis. Note that multiple knots can be used in order to induce discontinuities in the basis functions.
Reproduced from [53].

references are discussed in the review. To help the reader we have divided the papers into
sections although some papers in practice cover more than one category. The general references
havebeen given separately for each section to reduce the time needed to find a particular
For strong-field problems, the B-Spline basis is a popular choice. B-Splines, due to theirreference. If a reference is mentioned in several chapters it is included only once (the firsttime). calculable matrix elements, are efficient in both memory and computation. The BanalyticallyThroughout, we use atomic units (au) unless specifically indicated.
Spline basis functions are used as is (as shown in Equation (3.20)), or as a basis to represent the

2. Description of B -splines

field-free eigenstates of the system.

B-splines are functions designed to generalize polynomials for the purpose of approximating
arbitrary functions. One can view them as new elementary functions such as sin x or jl (x).
One becomes familiar with them by understanding their qualitative behaviour and how to usethem, that is how to obtain values of the functions, their derivatives or integrals. This is the

36
TDSE calculations of the ATI spectra of hydrogen [55] was one of the first successful application of B-Splines to the strong-field regime. Since then, B Splines have been used to examine highharmonic generation, fully two dimensional multi-photon ionization of helium, and examinationsof single ionization in H− . For more information on the applications of B-Splines to atomic physics
calculations, see [53, 56].
Slater type orbitals are another popular basis set for radial calculations in atomic and molecular physics [47, 57]. They are typically of the form:
Sn,ζ (r) = rn−1 e−ζr ,(3.26)

where n and ζ are parameters of the basis function. Typically, n is associated with the principal
quantum number, while ζ is a tunable parameter which describes the radial extent of the orbit.
Slater orbitals were originally investigated due to their similarity to the hydrogenic wavefunctions.
The Slater orbitals form a complete basis for a fixed ζ. However, slow convergence can occur forfeatures much farther from the origin than ζr. In order to combat this, it is common for multiple
values of ζ to be used in a single calculation, expanding the size of the basis.
Slater type orbitals are typically used in variational calculations of atomic states [57, 58].
This allows basis sizes to be minimized while still maintaining accuracy for single (or few) statecalculations. Unfortunately, this variational process produces a new set of basis functions for eachnew state, which makes calculations of large numbers of states difficult. As the number of nodes
in a wavefunction goes up, the number of basis states required also increases. The combination ofthese two behaviors – the individualized basis states for each state, and the increasing number ofbasis states necessary for highly oscillatory states – means that Slater type orbitals have not seenmuch use in strong-field calculations.
Gaussians are another common basis set used in solutions of the atomic and molecular
Schrödinger equation [59]. There are a number of different schemes for using Gaussian functions
as a basis for finding atomic orbitals [47]. A common one is analogous to the Slater type orbitals:2Gl,α (x) = rl e−αx ,(3.27)

37
where α is a parameter defining the width, and l defines the order of the basis element. For a fixed
α, the set of Gaussians of different order l forms a complete basis. Unfortunately, a basis set with
fixed α is exceedingly slow to converge for most realistic atomic orbitals, and thus multiple basis
sets with different α’s are used to speed convergence.
Because of this, calculations using Gaussian basis sets are usually made up of composite
basis elements, where a linear combination of pure Gaussian basis elements create a set of new
basis elements. A common choice is using linear combinations of Gaussians to generate Slater typeorbitals. Because Gaussians are simpler to analytically integrate than Slater type orbitals, thismakes calculations with Slater type orbitals easier to manage, at the expense of accuracy, as afinite number of Gaussians cannot represent a Slater orbital exactly.
Slater-type orbitals and Gaussian-type orbitals are commonly used in computational chemistry. They can generate very quick, accurate electronic structure calculations of multi-centeredmolecules due to their natural similarity to real atomic orbitals. Unfortunately, both suffer fromdifficulty in accurately representing the continuum wavefunctions. This inability to describe thecontinuum efficiently has decreased their usefulness in strong-field calculations, though some hybrid calculations have been able to utilize these wavefunctions for bound state representations,combining them with more flexible methods for continuum state calculations [60].

3.2.2

Numerical basis state method

We find the field-free basis numerically by solving the eigensystem given by Equation (3.12).
The spectra of any unrestricted atom is infinite in the number of bound states as well as the numberof continuum states. Of course, this cannot be represented in finite memory, nor calculated withfinite computational power. In order to overcome this, we restrict the atomic spectra by confining
the atom within an infinite spherical well potential, which discretizes the continuum part of theatomic basis states, and reduces the number of bound states to a finite number. We artificiallyrestrict the number of continuum states by restricting the maximum energy state that we need torepresent in our wavefunction.

38
In practice, the range of the infinite spherical potential is denoted by Rmax , and a Dirichletboundary condition is enforced at the boundary, i.e. Φ(Rmax ) = 0. The value of Rmax sets thedensity of states for the continuum and highly excited bound states. We typically use values of
500 − 1000 a.u. for our calculations (see section 3.2.3 for more information on how this affects ourbasis set.)
As noted above, we use the spherical harmonics as angular eigenfunctions, while the radial
eigensystem is given by Equation (3.11). We use a finite difference method to find the individual
eigenfunctions represented as values on a radial grid. In order to maintain accuracy around thesingularity in Equation (3.11), a higher density of points is needed near the singularity, while
far fewer are required far from the singularity. Unfortunately, finite difference calculations withnonlinear grid spacing are often plagued by numerical errors. Transforming Equation (3.11) intoa logarithmic variable (denoted by x) gives us the required accuracy at small values of r andincreases the spacing at higher values of r, while allowing us to maintain an evenly spaced grid.
Such a nonlinear grid is accomplished by transforming Equation (3.11) using x = log(r), and wepfind the transformed wavefunction ηn,l (x) = Rn,l (r(x)) r(x):#"∂xx (l + 12 )2
2(3.28)++ r(x) (En,l − V (r(x))) ηn,l (x) = 0.
−22Here, the x coordinate is evenly spaced, while the r coordinate is logarithmically spaced. We
have removed the explicit notation of our boundary conditions from the potential, these will be
enforced by the numerical method rather than as a term in the potential. This necessitates theneed to set a minimum radial distance (denoted Rmin ). In test calculations we have found that
Rmin = 10−6 a.u. leads to accurate representations of the ground state wavefunctions. We constrain
this second order differential equation to have a unique solution by fixing the boundary conditions
Rn,l (Rmin ) = Rn,l (Rmax ) = 0, and ensuring the resulting function Rn,l (r) is normalized.

39
3.2.3

Properties of the basis

There are a number of free parameters that we must decide on when calculating our basis.
Here we list the parameters, what impact they have on the calculations, and typical values that wehave used in our calculations.
• Rmax : the largest r value in the grid, defined by the boundary conditions. This determines
the spatial extent of the wavefunction and the density of the continuum states. Figure 3.2
shows the energy vs. n of the state for different Rmax . The quadratic dependence ofthe energy as a function of n is fairly clear, along with the increased density of states(decreased quadratic coefficient) as Rmax increases for fixed n. Unlike typical grid based
methods, we do not have an absorber at Rmax , so we require a very large grid, to avoidreflection problems. The strength of the electric field we can represent is limited by thishard boundary: any part of the wavefunction that reaches the boundary of our grid will bereflected off. Test calculations have shown that a good lower bound on the extent of thegrid with respect to electric-field strength E and frequency ω is:
Rmax & N

πp
20UP ,
ω(3.29)

where N is the number of periods of the laser field, ω is the frequency of the laser fieldand Up is the ponderomotive potential E 2 /4ω 2 . An electron that is ionized near the centerof the pulse and excited to 10Up (the maximum classical rescattering energy possible in asingle cycle [61, 62]) will have just enough momentum to drift to the edge of our box by theend of the pulse. For a 1014 W/cm2 and 800 nm laser, this corresponds to approximately
130 a.u. per 2 cycles of the laser field. We typically use 1000 atomic units as our maximum
radius.
• npoints : The number of points on the grid. Due to the logarithmic spacing of our grid
points, the density is lowest at Rmax . This restricts the maximum n state that we canrepresent: higher n states have higher frequencies, and the Nyquist frequency requires that

40
Rmax
300
400
500
600
700
800
900
1000

120

Energy (a.u.)

100
80
60
40
20
0
0

200

400

600

800
n

1000

1200

1400

Figure 3.2: Energy as a function of n for different Rmax from 300 to 1000. Note the quadratic
dependence of En as a function of n that becomes more pronounced for smaller Rmax .

we have at least two points on our grid per period of the wave. Practically speaking, wehave found that significantly more points than that must be used to prevent errors both incalculating the states and in integrating the states to find the dipole transitions. We show
the effects of using too few points for the wavefunction in Figure 3.3, where the left panel
shows the high r part of the wavefunction when too few points are included, and the rightshows the same part of the wavefunction when we include enough points. The number ofpoints is functionally limited by the amount of memory available: we use long doubles (16bytes) for each point in our calculations, so we require at least npoints ∗ nmax ∗ 16 bytes ofmemory per node. Typically, we use 300, 000 points for nmax = 2000 and Rmax = 1000.
• nmax : the maximum principal quantum number state that we can represent, which is anupper bound on the energy represented by a given wavefunction. This parameter is limitedby the number of points on the grid and Rmax . A lower bound on the required nmax for an
atom interacting with an electric field E with frequency ω is 10Up = 10E 2 /4ω 2 . For a 1014
W/cm2 and 800nm laser, this corresponds to approximately 3 a.u. For nmax = 2000 and
Rmax = 1000, En=2000 ≈ 19 a.u. The relationship between Rmax , n and E is shown above
in Figure 3.2.

0.04

0.04

0.02

0.02
calculated

0.00

(r)

(r)

41

0.00

0.02

0.02

0.04

0.04

900

920

940
960
r (atomic units)

980

1000

900

920

940
960
r (atomic units)

980

1000

Figure 3.3: A demonstration of aliasing shown for state n = 2000, l = 0 with Rmax = 1000, with
100, 000 points on the left as an example of too few points, and 300, 000 points on the right as an
example of barely enough. For clarity, only the last 100 a.u. of the wavefunction has been shown.
Note the oscillating amplitude for the plot with too few points. This is an example of a state that
does not have enough points for an accurate integration when finding the dipole moment.

• lmax : the maximum angular quantum number we can represent. Due to the fact that weput our potential inside an infinite spherical well potential, we do not have the accidental
degeneracy of the l states [63], and our highest l states have more energy than the corresponding low l states with the same principal quantum number. Typically, we set lmax to70.

3.2.3.1

Representation in the field-free basis set

We are interested in representing the solution of the Schrödinger equation in the field-freeenergy basis. Using the radial eigenfunctions and spherical harmonics, any wavefunction can be
written as:Ψ(r, t) =X

cn,l (t)Rn,l (r)Yl,m (r̂),(3.30)

n,lwhere Yl,m (r̂) are the well known spherical harmonics found in e.g. [64], and m is constant for
any given calculation (typically m = 0). Rn,l (r) are the radial wavefunctions that we found in theprevious sections, and cn,l (t) are the time dependent complex coefficients. In the field free basis,
we represent the wavefunction Ψ(t) as a vector composed of the complex coefficients cn,l (t). Thus,a wavefunction denoting the ground state of hydrogen is merely |c1,0 (t)| = 1 and cn6=1,l = 0. We

42represent this vector as C(t).
The time dependent Schrödinger equation for an atom interacting with a oscillating electricfield pointing in the z direction can be written down as:


H0
z
}|
{ H 0 (t) p2z }| {+ V (r) + zE(t) Ψ(r, t) = i∂t Ψ(r, t).
−
 2


(3.31)

The field-free Hamiltonian H0 is diagonal in its own basis, and thus the projection of this term
into the field-free basis is trivial. The time dependent part of the Hamiltonian H 0 (t) can be foundwithin the field free basis as:ẑ

n,l,n0 ,l0= hφn,l | r cos θ |φn0 ,l0i=

Z RmaxRn,l rR0n0 ,l0drZ

Yl,m cos(θ)Yl0 ,m dr̂,(3.32)

where ẑ is the dipole moment operator in the field-free basis, with matrix elements given by the
integrals in Equation (3.32). In practice, the spherical harmonic integrals are found from theClebsch-Gordon coefficients using a library function [65] for more efficient calculation. Combining
this information, the time dependent Schrödinger equation is given by:hi
Ĥ(t)C(t) = Ĥ0 + E(t)ẑ C(t) = i∂t C(t),(3.33)

where we have defined the time dependent Hamiltonian in the field free basis as Ĥ(t).

3.2.4
3.2.4.1

Implementation of numerical basis state method
Numerov Method

We would like to solve Equation (3.28) numerically on an evenly spaced grid, with respect tox. However in this form, it is no longer an eigensystem of the form:Âv(x) = λv(x),(3.34)

where Â is a matrix, and λ and v(x) are the eigenvalue and eigenvector respectively. Instead,
Equation (3.28) is known as a generalized eigenvalue problem, which are equations of the form:Âv(x) = λB̂v(x),(3.35)

43
where B̂ is also a matrix. Unfortunately, in our case, B̂ = r2 is a matrix that is very poorly
2 /R2 . Typical values of κ for ourconditioned, with an approximate condition number of κ = Rmax
minproblems are 1018 , since the loss of precision in basic arithmetic methods is typically the log of thecondition number, log10 (κ) = 18. Finite-difference matrix based eigensolvers become unrealistic toapply to this formulation of the problem.
Instead, we use the Numerov method, which is an explicit, fourth order accurate, linear
multistep method for solving linear second order differential equations. The Numerov method canbe applied to problems of the form:∂xx y(x) = f (x)y(x),(3.36)

where y(x) is the solution to the second order differential equation to be found, and f (x) is thedriving function. Equation (3.28) can be written in this form with:f (x; En,l ) =(l + 21 )2+ r(x)2 (En,l − V (r(x))) ,2(3.37)

where it is important to note the dependence of f (x) on the energy of the state En,l .
The Numerov method is part of a family of methods known as “shooting” methods. Thesemethods use a backward difference to find the next value of the solution y(x) from the previous twovalues already found. In the case of the Numerov method the next value of the solution is given
by:





2
h2
y(xn ) 1 − 5h
f
(x
)
−
y(x
)
1
+
f
(x
)
n
n−1
n−1
12
12


y(xn+1 ) =,5h2
1 + 12 f (xn+1 )(3.38)

where xn is the nth point in the continuously spaced grid between Xmin = log(Rmin ) and Xmax =log(Rmax ).

3.2.4.2

Algorithm

It is important to note that the Numerov method is not an eigenvalue solver but a differential
equation solver. It attempts to find a function y(x), given a function f (x; En,l ) and two starting
points. However, for an energy that is not an eigenvalue, y(x) will eventually diverge. We must use

44this divergent behavior – and the other properties of the eigenfunctions such as continuous second
derivatives for continuous potentials and adherence to the boundary conditions – to determine analgorithm that will reliably give us accurate eigenvalues and eigenfunctions.
The algorithm that we use is fundamentally a binary search. We first establish an upperbound and a lower bound to the energy of the eigenstate that we are interested in finding. We thencalculate the wavefunction for the energy in the middle of this region and determine if that energy
is greater or smaller than the energy of the eigenstate that we are interested in. If it is greater, we
use it as a new upper bound on the energy, and look for an energy halfway between this new upperbound and the previously established lower bound. If it is lower, we similarly restrict the search
space. In this way, we converge to the correct energy to whatever accuracy we desire, within thediscretization error of the Numerov method.
The ability to do this binary search relies on the ability to accurately determine whether theguess for the energy of the state is high in energy or low in energy in comparison to the true energyof the state. There are two different techniques available to determine if the state is high or low
in energy. However, they suffer from robustness issues depending on the energy, and how far theenergy is from the correct energy for the state that we are trying to find:• The number of nodes (zeroes) of the wavefunction. We know that solutions to the sphericaltime independent Schrödinger equation with quantum numbers n and l have n − l − 1
nodes. However, for tightly bound states, where the wavefunction approaches zero before
our external boundary condition must be enforced, there is a wide range of energies thatwill have the correct number of nodes. For energies below the ground state of the atom weare examining, the Numerov method can also exhibit behavior where it generates spurious
nodes.
• Enforcing the boundary conditions by propagating from both ends, and then determining
the derivative of the wavefunction where the two propagations meet. For an eigenfunction,these two wavefunctions should meet exactly, and have a matched first and second derivative

45
where they meet. However for an energy that is too small, the forward derivative will havea smaller second derivative than the backward propagation. For an energy that is too large,
the opposite occurs.
• When the solver is very close to the correct energy, the difference in the second derivative
can get smaller than machine precision, while there might be a difference between the upperbound and lower bound of the energy that is significantly larger. This necessitates the useof the first derivative to finish the calculation. For an energy that is too small, the forward
first derivative will be larger than the backward first derivative.
While none of these techniques can find a state on their own, we can combine them inorder to design a more robust algorithm for determining the correct eigenstate, with the correct
eigenenergy. After we calculate the wavefunctions, we normalize them, and store them for later
orthogonalization.
This processes is independent for each value of l that we wish to use, allowing for a simple
parallelization in calculating the eigenstates for each l independently. This results in each processor
working on a single l value, and solving for each n in sequence. Parallelization of the individual neigenstates is not worth the effort due to the sequential nature of the orthogonalization that mustbe performed after finding each n.

3.2.4.3

Orthogonalization

Solving for basis states using the Numerov method is not exact. Small numerical inaccuracies and instabilities prevent the wavefunctions determined this way from being pure eigenstates.
Furthermore, since we are not using an eigenstate solver that has an orthogonalization step, theorthogonalization must be enforced explicitly. This is important, since the dipole transition operator within the field-free Hamiltonian needs to be found from orthogonal states in order to avoid
errors in the calculation.

46
The states that are found by our eigensolver are a mixture of the true orthogonal eigenstates
φen,l (r) = an,l φn,l +X

n0 6=n,lhφen,l |φn0 ,l i φn0 ,l (r).

(3.39)

Here, φe represents a state from the eigensolver that still needs to be orthogonalized, φ represents apure state, and an,l is the fraction of φe that is a pure state. For the states found by our eigensolver,an,l will be close to 1.

We remove the unwanted contributions from each state using the modified Gram-Schmidt orthogonalization routine [66]. The modified Gram-Schmidt orthogonalization routine just subtracts
the projection of the completed states from each new state after it is calculated. In this manner,
only the components of the previously calculated pure states are left:
Algorithm 1 The Modified Gram-Schmidt algorithm
1: for all n ← nmin . . . nmax do
2:
vn ← φen,l3:
for all i ← 1 . . . n − 1 do
4:
vn = vn − (vi · vn )vi5:
vn = vn /|vn |Algorithm 1 is only concerned with a growing subset of the basis states in our basis set –the outer loop (line 1) is done concurrently with the Numerov eigensolve. Because each state is
only orthogonalized against the previous states that have been calculated, we can orthogonalize astate in parallel while finding the next state via the Numerov method. This significantly reduces
the calculation time.
Note that the orthogonality condition only needs to be enforced for the radial eigenfunctions
with the same l value. States that have different l values have their orthogonality ensured by the
angular part of the wavefunction.

3.2.4.4

Comparisons with previous results and literature data

Any numerical method should be compared against known results. Thankfully, the hydrogen
atom is analytically solvable, and gives us a good benchmark to compare against.

47
The analytic solutions of the radial states of the hydrogen atom can be written down in termsof Kummer’s function:
iRn,l (ρ) = N −
ρ
2Zr
ρ=
,
n

l

e−ρ/2 M (l + 1 − n, 2l + 2, ρ),(3.40)(3.41)

where N is a normalization constant, Z is the nuclear charge, n is the principal quantum number
for bound states, and n = −iZ/k for continuum states with momentum k.
For the bound states, this reduces to:
Rn,l (r) = N e

−r/n



2r
n

l+1

L2l+1
n−l−1



2r
n



,

(3.42)

where L is one of the associated Laguerre polynomials, and each state has energy En = 1/(2n2 ).
The continuum states are given in terms of the Coulomb wave function (Fl (η, ρ)):
−1Rn,l (r) = N Fl
, kr ,
k(3.43)

which we find using a library function [67, 68]. The momentum is calculated from the dispersion√
relation k = 2E.
A comparison between our calculations and the results of Equations (3.42) and (3.43) has
been performed. In Figure 3.4, the bound states are plotted along with the relative error between
the analytic and calculated wavefunctions, while in Figure 3.5, the same has been done for thecontinuum. The boundary conditions change the character of the bound states with principal
quantum number above roughly n = 16, as the exponential decay can no longer get small enoughbefore encountering the infinite potential at the edge.
For the bound states, we can also compare the calculated energies to the analytically known
energies. In Figure 3.6, the relative error is shown for the first five l eigenvalues and the first 18
n eigenvalues. The calculated eigenvalues diverge from the real eigenvalues drastically when n isgreater than 16, as the boundary conditions influence the wavefunction and its energy.
We have adapted this method from a previous work by Chen et al. [69], but we improved
upon the accuracy of the bound states found in [69] significantly. Table 3.2 gives a comparison

48
calculated
analytic
(r)

0.4

0.0

0.2

0.1

0.0

0.2

10 12
10 15
10 18

0

2

4
6
r (atomic units)

8

calculated
analytic

0.1

rel. error

rel. error

(r)

0.6

10

10 6
10 10
10 14
10 18

0

25

50

75
100
125
r (atomic units)

150

175

200

Figure 3.4: A comparison between the analytic and the calculated bound states for the ground state(n = 1, l = 0) on the left, and n = 5,l = 1 on the right. The relative error of |Ψcalc − Ψanal |/Ψanal isincluded underneath each plot. These were calculated with Rmax = 1000, 300, 000 points per grid,and Rmin = 10−6 . The plots only show a portion of the r range of the wavefunctions so details canbe seen.

between this work and the previous work. We note that we have reached an improvement in therelative error between the analytic result and our calculated energies on the order of 10−3 to 10−4
better.

3.3

Time dependent calculations
The ab-initio calculation of the susceptibilities requires an ab-initio calculation of the timedependent wavefunction. Thus, we must propagate Equation (3.33) in time, which usually takes
the vast majority of the calculation time. There are two important parts of the numerical solutionof the time-dependent calculations: the propagation itself, which is done using the Crank-Nicholson
method for the work presented in this thesis (see section 3.3.1) and the implementation of an absorbing boundary condition, which in this work is related to the absorption of part of the wavefunction
representing high energy states.

3.3.1

Propagation calculations

Equation (3.33) can be approximated as:
C(t + ∆t) u e−iĤ(t)∆t C(t) u

1 − i ∆t
2 Ĥ(t)1 + i ∆t
2 Ĥ(t + ∆t)C(t),(3.44)

49
0.04
calculated
analytic

0.00
0.02

(r)

(r)

0.02

0.00
0.05

10 7

rel. error

rel. error

0.04

10 12
10 17

calculated
analytic

0.05

0

200

400
600
r (atomic units)

800

1000

10 3
10 7
10 11
10 15

0

200

400
600
r (atomic units)

800

1000

Figure 3.5: A comparison between the analytic and the calculated continuum states for the staten = 30,l = 1 on the left, and n = 200,l = 70 on the right. The relative error of |Ψcalc − Ψanal |/Ψanalis included underneath each plot. These were calculated with Rmax = 1000, 300, 000 points per
grid, and Rmin = 10−6 . The momentum was found from the calculated energy of the state, rather
than the quantum number. This compounds the error in the wavefunction, and results in a large
error near the boundaries as the analytic wavefunction is not constrained to be zero due to theboundary conditions.

for a small timestep ∆t. The second step uses the Cayley unitary transformation [70] in orderto preserve the normalization of the wavefunction. The representation of the time propagation
operator via the exponential function leads to the following condition: the exponent Ĥ(t)∆t has tobe much less than 2π, in order to guarantee that we do not get aliasing effects from the propagationof the wavefunction. We satisfy this condition when Enmax ∆t < 1. This typically results in time
steps that are approximately 0.05 a.u.
This approach is known as the Crank-Nicholson method which is a unitary, energy conserving,
and unconditionally stable method for propagating first order in time differential equations tosecond order accuracy [71]. For our system, we rewrite Equation (3.44) as:





∆t
1 + i ∆t
Ĥ(t
+
∆t)
C(t
+
∆t)
=
1
−
i
Ĥ(t)
C(t).
2
2

(3.45)

∆tThis is an equation of the form Ax = b, with A = 1 + i ∆t
H(t
+
∆t)
,
b
=
1
−
i
H(t)
C(t)
2
2

and x = C(t + ∆t).

The algorithms just described can be formulated in terms of repeated solutions of linear
matrix equations. Due to the ubiquity of this method, a number of well designed software packages

50

10−13

10−1

10−14

10−3

10−15

10−5

10−16

10−7
0

5

10

15

10−9

rel. error

rel. error

10−12

1
0
1
2
3
4
0

5

10

15
n

20

10−11
10−13
10−15

25

Figure 3.6: The relative error in the calculated energy eigenvalues compared to the analyticalenergy En = 1/(2n2 ) for l < 5. The right figure is a zoomed in version of the left, to better show
the errors at low n. The n > 16 states start to be influenced by the bounding box, while n > 23
states have errors larger than 10%, shown by the dot-dashed brown line.

are available to help in the solution. We solve this linear system using PETSc: the Portable
Extensible Toolkit for Scientific computation [72–74].

3.3.2

Absorbing boundary

Strong field ionization processes can generate very high energy electrons upon the absorptionof a large number of photons. Since we need to restrict our calculation to a finite number of basis
states, we can only represent states up to some maximum energy. Consequently, we need to suppress
the population of even higher energy states. We solve this problem with an absorbing boundary
implemented for states with either high energy (high n states) or high angular momentum (high lstates). The absorbing boundary is represented by a masking function which is multiplied to the
wavefunction at every time step. The masking function used in this work is a cos1/8 function givenby:
1/8

a(n, l) = cosn − (nmax − nabs ) π
nabs
2



cos

1/8



l − (lmax − labs ) π
labs
(2)



(3.46)The specific function used in the present work was found to sufficiently reduce the wavefunction
inside the absorbing region and, hence, prevent reflections back to states with lower energy and/or
angular momentum. In test calculations we have found that in this specific set of problems thecos1/8 absorber leads to a more efficient suppression of the wavefunction than other techniques

n
1
2
5
10

En
−0.5
−0.125−0.02−0.005Prev. Result
l=0
l=1
−8
2.0 × 102.0 × 10−8 0.0 × 10−10
8.0 × 10−9 5.0 × 10−9
2.3 × 10−9 1.9 × 10−9This Work
l=0
l=1
−14
4.3 × 10
2.6 × 10−13 4.6 × 10−14
2.1 × 10−14 3.9 × 10−14
9.2 × 10−14 3.3 × 10−15

51

Table 3.2: The relative error between the analytic and calculated energies of a few different boundstates from [69] and this work. We achieve significantly better (3 or 4 orders of magnitude) accuracy
using our logarithmic grid.

such as complex absorbing potentials or other potentials of the form cos1/n . In Figure 3.7, the final
wavefunctions of a 10 cycle sine-squared pulse with carrier frequency 800 nm and intensity 1014
W/cm2 with and without an absorber.

3.4

Summary
In summary, the numerical basis state method is a capable process for doing ab-initio calcu-

lations of the time dependent Schrödinger equation. The field-free basis states that underlie boththe numerical basis state method and the perturbative calculations that will be used in subsequent
chapters are found using a Numerov based eigensolving algorithm, while the propagation used in
the ab-initio is performed using the Crank-Nicolson method.
This will allow us to compare between the perturbative and non-perturbative calculations of
the susceptibility, as seen in Chapter 4. Furthermore, we will use this method on its own in order
to find the time dependent susceptibility in Chapter 5.

52

10−6

10−7

−10

10−11

10−14
10−18

10−15
10−19

−22

10−23

10−26

10−27

10−30

10−31

10

50 n states, 20 l states
200 n states, 20 l states
500 n states, 20 l states

10−3

population

population

10

101

50 n states, 5 l states
50 n states, 10 l states
50 n states, 20 l states

10−2

0.0

2.5

5.0

7.5

10.0
energy (a.u.)

12.5

15.0

17.5

0.0

2.5

5.0

7.5

10.0
energy (a.u.)

12.5

15.0

17.5

Figure 3.7: The final wavefunctions from two different calculations using a 10 cycle sine-squared
pulse with carrier frequency 800 nm and intensity 1014 W/cm2 . Panel on the left shows threedifferent sizes of l absorbers for an n absorber of size 50. Good convergence is achieved, with l = 20
giving good accuracy at lower energy levels. Panel on the right shows three different sizes of nabsorbers with an l absorber of size 20. All absorbers are identical for lower energies, however then absorber of size 50 shows some reflection at high energies.

Chapter 4

Calculations of ab-initio and perturbative susceptibilities

In this Chapter we apply the numerical basis state method to calculate the nonlinear perturbative susceptibilities, and compare them to ab-initio calculations of the susceptibilities. Throughthis comparison, we show that the perturbative susceptibilities we calculate fail to show converging
behavior at intensities above a certain threshold. We further find that in the same intensity regime
the results of perturbative calculations and the ab-initio calculations diverge. We examine a few
different systems in order to confirm this behavior, and discuss the effect that the ionization energy
has on divergence as a function of intensity using hydrogen, helium, as well as different model
atoms.
We also examine the process of low-order harmonic generation by examining the third andfifth harmonics in hydrogen. By analyzing the perturbative series as a function of intensity, andthe ab-initio calculations as well, we determine the breakdown of the perturbative approximation
for low order harmonic generation.

4.1

Introduction
The nonlinear response of an atom to a strong external field, e.g. an intense laser, can giverise to a number of interesting phenomena. Low order harmonic generation, where the atom emits
odd harmonics of the fundamental frequency, is generated by absorbing n photons at frequency ωand emitting a photon with frequency nω. If nω is large enough to ionize the atom, multiphoton
ionization is observed. These processes are perturbative at low intensities: the efficiency of nω

54
generation is proportional to E n , where E is the electric field strength.
The response of the polarization to the electric field is given by:
P (ω) = χ(ω)E(ω),(4.1)

where χ is the electric susceptibility of the material. For a perturbative process, this can be
expanded in powers of the electric field as:
P (ω) = χ(1) (ω)E(ω) + χ(2) (ω; ω1 , ω2 )E(ω1 )E(ω2 ) + χ(3) (ω; ω1 , ω2 , ω3 )E(ω1 )E(ω2 )E(ω3 ) + . . . (4.2)
where ω =Pi ωi .

Typically, this series is truncated after the first nonzero nonlinear term χ(3)(χ(2) is zero for centrosymmetric media), and the term where ω = ω1 = ω2 = −ω3 can be used todescribe the physics of self-focusing.
We can see the three step process of filamentation (discussed in Chapter 1) demonstrated in
the ab-initio calculations of the susceptibility of noble gases done in 2002 by Nurhuda et al. [75].
In Figure 4.1 we show the susceptibility as a function of intensity. The linear increase in thesusceptibility at low intensity is self-focusing, the saturation is associated with the ionization, andthe downturn is associated with plasma induced defocusing. This saturation and eventual downturn
requires negative higher order contributions in the susceptibility, though Nurhuda explained thesaturation and downturn as the influence of the ionized electron.
The ionization based model described in Chapter 1 was believed to be well understood [15].
However, in 2009, Loriot et al. challenged this model by extracting out higher order terms of theindex of refraction [24] (see Section 2.1 for how the index of refraction is related to the susceptibility)
by using a cross polarization measurement that is insensitive to the plasma contribution. The
extracted terms, part of the perturbative series for the index of refraction:n = n0 + n2 I + n4 I 2 + n6 I 3 + n8 I 4 + n10 I 5 + . . . ,(4.3)

shown in Table 1.1 have a positive n2 , a negative n4 , a positive n6 and a negative n8 . This challenged
the idea that the saturation and downturn of the susceptibility was caused by ionization alone, andled to the Béjot et al. paper where they examined filamentation without ionization [25].

55

Figure 4.1: Susceptibility of the neon atom as a function of intensity. Note the linear increase as afunction of intensity for low intensities, the saturation, and the eventual decrease of the suceptibility.
From [75].

This possibility subsequently ignited a significant controversy in the field of filamentation,and there have been a number of papers that have attempted to either confirm or disprove theapplicability of the higher order Kerr effect (HOKE) contribution to the filamentation process. It
raised two related questions: Are the higher order Kerr terms an accurate way to describe thephysics of filamentation? Did Loriot et al. actually extract out the perturbative Kerr terms? We

Fig.notable
1. theoretical
Plot of considerations
the real part
ofhigher
nonlinearsusceptibility∆χ owill discuss a few of the moreof theorder Kerr
terms.

by direct numerical integration of TDSE as a function
intensity. In panels a), b), c), and d), results for Ne,
the hydrogen atom [29]. They attempted to extract the susceptibilities in two different ways: by
and Xe are shown. The laser frequency used for simulati
examining the nth harmonic yieldas aeV.
proxy for χ(n) (the harmonic susceptibilities), which can be
1.55
Köhler et al. used an ab-initio method to determine the higher order susceptibilities ofdone as a function of intensity; and by finding the individual perturbative degenerate susceptibilities
from the susceptibility vs. intensity curve. For a perturbative process, the susceptibility terms

where it has been found a cusp in low intensity domainincrease in low intensity regime is in accordance with t
the nonlinear susceptibilities of hydrogen are reproduced from [29]. Notice that for longer pulses, theturbation theory prediction. It is also pointed out thatsusceptibilities are all constant, or nearly so at low intensities. This is the behavior that is expected
the linear increase of susceptibility, only transitions b
for perturbative susceptibilities. However, as the intensity increases, the susceptibilities that are
ground state and excited states are important for th
extracted are no longer independentof the intensity,
implying
a fundamentally
non-perturbative γ have r
cesses. Several
works
devotedfor calculatingthat it is suﬃcient to include the bound-excited statsitions,17) thus neglecting the role of bound-continuum
sitions. In connection with the nonlinear susceptibili
should be constant with intensity by definition. However, this is not what was found. In Figure 4.2,

process.

approach evaluates the susceptibilities at different harmonicfrequencies j ω0 , whereas the latter approach gives access toall χ (j ) directly at ω0 . For purely instantaneous HOKE termsof Eq. (5), both approaches should yield the same χ (j ) .
Let us first evaluate χ (j ) at the harmonic frequencies
j ω0 . Results are presented in Figs. 1(a)–1(c), where χ (3) ,(5)
χ , and χ (7) are plotted for different pulse durations τ vspeak intensity. We can observe how difficult it is to extract
meaningful coefficients χ (j ) being independent of intensityand pulse duration in the range where PNL becomes relevant,i.e., I > 1 TW/cm2 . The fact that it is not possible to describe
the nonlinear polarization in terms of a Taylor expansion at
intensities above 10 TW/cm2 already indicates significant
changes in the electronic configuration of the atomic system.

strongly fluctuate along the intensity range. We will further
refer to this simple model as HOKE_35. Note that we will
later introduce another set of parameters, HOKE_FIT, which
includes nonlinear coefficients up to χ (9) , chosen to obtain amore accurate description of the nonlinear saturation.
For the standard model with MPI rate, in addition toχ (3) , the photon number K and the cross section σK has tobe determined from the TDSE simulation data. For pulses
with wavelength λ0 = 800 nm and the ionization potentialof hydrogen of UH = 13.6 eV one obtains a minimum photon
number of K = mod(UH /h̄ω0 + 1) = 9. However, we find theleast dependence of σK on the peak intensity when employing
a reduced photon number of K = 6. This is plausible, sincefor the intensities employed here the Keldysh parameter γ
−45

x 10

x 10

8 (a)

0 (b)

4

((m/V) )

10 fs
20 fs
100 fs

6

−5

(5)

4

χ

2

χ(3) ((m/V) )

−25

56

2
−1

0

10

10

−10 −110

1I (TW/cm2)

10

x 10
(d)

6

6

4

2

σ6 ((cm /TW) /s)((m/V)6)
(7)

1

10

I (TW/cm2)

−30

x 10
8 (c)

−65

χ

0

10

2

1

0.5

0
−2 −1
10

0

10

0
0

1

2

10

I (TW/cm )

20

40
2
I (TW/cm )60

FIG. 1. (Color online) Model parameters extracted from 1D TDSE simulations. (a) χ (3) , (b) χ (5) , (c) χ (7) for τ = 10 fs (blue circles),from 1Dτ =nonlinear
20 fs (red squares) andτ = 100 fs (green
crosses) as a function of theof
peaktheintensity
I . (d) The MPI crosssection as
σ extracted
Figure 4.2: The
harmonic
susceptibilitieshydrogenatom
a functionof intensity.
TDSE simulations [same color coding as in panels (a)–(c)].
From [29].
6043811-3

Bree et al. attempted to calculate the nonlinear indices in a different way [31, 32]. They first
calculated the n photon ionization cross-section using a modification to the strong-field approximation, and used this as a proxy for the nth order perturbative absorption. Using the nonlinear
Kramers-Kronig relationship, they found the approximate nonlinear coefficients of the index ofrefraction. This index of refraction qualitatively showed the same saturation and downturn of thesusceptibility present in the Loriot et al. paper, however, theoretical questions about the use of aionization cross section to find an ionization independent index of refraction lay doubts as to thevalidity of this result.
These theoretical and experimental studies hinted at, but failed to reach, a consensus on
the applicability of the higher order susceptibility terms to the physics of high intensity pulses.
We intend to determine the boundary of applicability for the perturbative susceptibility with ourwork, by making a comparison between perturbative and ab-initio calculations within the same
theoretical framework.
The remainder of this Chapter is spent discussing the results of the comparison between thecalculated perturbative susceptibility and the calculated ab-initio susceptibility (section 4.2) and57
an examination of the perturbative and ab-initio low order harmonics (section 4.4).

4.2

Electric Susceptibility
We first examine the electrical susceptibility. We would like to discuss the higher order Kerr

effect, however, the index of refraction is fundamentally a macroscopic phenomena. We discussed
in section 2.1 the relationship between the index of refraction and the electrical susceptibility1
Any discussion of the validity or accuracy of the perturbation series requires a benchmark tocompare it to. Thus, our discussion of an intensity limit for the higher order perturbation terms
requires a comparison between the perturbative series and ab-initio calculations. Furthermore, inorder to reduce the relevance of various numerical approximations, we use the same basis states tocalculate both the ab-initio and perturbative calculations.
We have discussed the calculation of the basis states necessary for both the perturbative series
calculation and the ab-initio calculation in subsection 3.2.4. We have also previously discussed theab-initio propagation in the field-free basis in subsection 3.3.1.

4.2.1

Ab-initio calculations

From Equation (2.12), we find the susceptibility by first determining the expectation valueof the dipole moment as a function of time for an electric field polarized in the z direction:
hp(t)i = −e hΨ(t)| ẑ |Ψ(t)i ,(4.4)

where |Ψ(t)i is the computed ab-initio wavefunction as a function of time, and ẑ is the dipole
moment operator given by:ẑ =X

n,l,n0 ,l0hφn,l | z |φn0 ,l0 i ,(4.5)

where |φn,l i is the field free energy state with principal quantum number n and angular quantum
number l.
1

Part of the results presented in this section are published in A. Spott et al., Ab-initio and perturbative calculationsof the electric susceptibility of atomic hydrogen, Physical Review A 90, 013426 (2014).

58
The electric field is represented as a sin squared envelope:E(t) = E0 sin2



πt
T0



sin(ωt)(4.6)

where E0 is the maximum field strength, T0 is the duration of the pulse, and ω is the central field
frequency. The susceptibility is then calculated using the relation:
χ(ω) =F [hΨ(t)| z |Ψ(t)i]
E(ω)

(4.7)

Ab-initio calculations are typically done with a nmax of 2000, and a lmax of 70, which takes
approximately 8 hours on 20 nodes of the Janus supercomputer. Each node has 12 cores.

4.2.2

Perturbative calculation

For the perturbative calculations, we use the eigenstates of the field-free energy basis to calculate the N th coefficient of the perturbative power series expansion of the ground state wavefunctionin the external field, as given in Chapter 2:|ψ (N ) (ω1 , . . . , ωN )i =
X

jN 6=j0

···

(4.8)


N

−iω t
X 
Y hψji | µ̂E(ωi )e i |ψji−1 i  |ψjN iiPj1 6=j0 i=1
ωji − ωj0 −ωkk=1

where ωj0 is the ground state energy, ωk and ωi are the participating frequencies of the electric field,
ji denotes the state in the numerical basis set, and µ̂ is given by Equation (4.5). The lifetimes of theexcited states are neglected since all calculations performed in this study are far from resonance.
The N th-order term in the expansion of the single atom polarization in an overall nω process canthen be written:

where n = 1, 3, 5, . . .,

P

D

P

(N )

N
E
X
0
0
(nω) = Phψ (j ) | µ̂ |ψ (N −j ) i ,(4.9)

j 0 =0j ωj = nω and ωj = ±ω. P refers to the average of all permutations of thefrequencies. The symmetry of the electric field with respect to positive and negative frequency components (due to the fact that the electric field is a real quantity) allows us to rewrite Equation (4.9)

59
as:
D

(N )

P

(N )

N
E
Y
(N )
(nω) = 0 χnωE(ωi )(4.10)

i=1

with χnω is the N th-order term of the susceptibility at frequency nω due to contributing electricfield frequencies ωi . This allows us to calculate the susceptibility for any given harmonic n, frequency ω and order N . For this section, we are concerned with situations with n = 1, specificallythe degenerate N wave mixing term for a perturbative susceptibility of order N .

4.2.3

Results of perturbative series and ab-initio calculations for electric susceptibilities

Here, we first present the results of our calculations for the perturbative coefficients in the
power series expansion of the electric susceptibility showing the convergence of the results withrespect to the box size of the numerical grid R0 and the basis size. We then proceed by analyzing
criteria for the breakdown of the series expansion at high field strengths and identify a regime fora critical field strength. Finally, we compare the results of the perturbative calculations with thoseof the ab-initio calculations, and determine the regions of correspondence and divergence of these
two methods.

4.2.3.1

Perturbative higher-order nonlinear susceptibilities

In Figure 4.3(a) we present the results of our perturbative calculations for the susceptibility
coefficients of the power series for the hydrogen atom, Equation (4.2), up to χ(11) as a function of
the frequency of a uniform external electric field. Over the frequency regime shown the coefficientsare almost constant, but the higher-order coefficients tend to increase at the largest frequencies
considered in the present calculations. The increase indicates that one approaches the regime ofresonant transitions. When we examine the same figure for the helium atom, Figure 4.3(b), wenote that the resonant transition regime is at a higher frequency, as expected due to the larger
energy gap between the ground state and first excited state in helium. Since the life times of theexcited states are neglected in the present perturbative calculations, we have restricted ourselves

60

1010

χ(11)

108

χ(9)

106

ǫ0 χ(N ) /ρ (a.u.)ǫ0 χ(N ) /ρ (a.u.)

1014

χ(7)

106

χ(5)
χ(3)

102
0.00

104
102

χ(9)
χ(7)
χ(5)
χ(3)

χ(1)

0.01

χ(11)

0.02 0.03 0.04
frequency (a.u.)

0.05

0.06

(a)

100
0.00

χ(1)

0.01

0.02 0.03 0.04
frequency (a.u.)

0.05

0.06

(b)

Figure 4.3: a) Results of calculations for the different terms in the power series expansion for thehydrogen atom, Equation (4.2), as a function of frequency ω of a uniform electric field. b) Samefor the helium atom. Note the difference in the magnitude of the susceptibilities for hydrogen andhelium. The dashed line denotes 800 nm, the wavelength used for the majority of our calculations.

to calculations up to χ(11) and the frequency regime below 0.06 a.u. (corresponding to wavelengths
of approx. 800 nm and longer), in order to avoid the effects of the lifetimes on our calculations. We
may note here already that based on the perturbative results in Figure 4.3(a) and Figure 4.3(b) it
is unlikely that a negative slope in the nonlinear susceptibility, such as one seen in Loriot et al. [24],or the negative slope that Nurhuda shows in [75] can be explained via the power series expansion
in the present theory, since the first five coefficients of the series are found to be positive over thefrequency regime studied.
The strength of the higher order terms for the helium atom increase at a slower rate comparedto hydrogen. As we increase the order, the values of the susceptibility terms for helium are muchsmaller than for hydrogen. For example, χ(11) is approximately 107 a.u. for helium, while itis approximately 1013 a.u. for hydrogen. This can be explained by the larger energy difference
between the first excited state and the ground state for helium compared to hydrogen, which leads
to a weaker coupling between the states in helium.
The convergence of the results presented in Figure 4.3(a) was tested well with respect to the

1025

1022

1021
|ǫ0 χ(N ) (mω)/ρ| (a.u.)ǫ0 χ(N ) /ρ (a.u.)

61
1026

1018
1014

χ(11)

1010

χ(9)

10

χ(7)

6

χ(5)
χ(3)
χ(1)

102
10−2

0

100

200

300

400 500
nmax

(a)

600

700

1017
1013
109
10

N = 11
N =9
N =7
N =5
N =3
N =1

5

101
800

10−3

0

200

400

600 800 1000 1200 1400
nmax

(b)

Figure 4.4: a) Perturbative susceptibility coefficients χ(N ) in the power series expansion, Equation (4.2), as a function of the maximum principal quantum number nmax in the basis set for thehydrogen atom. The radial box size was R0 = 500. The dashed line denotes the transition frombound to continuum sets in the present basis set, while the crosses indicate convergence of theresults within 1% of the respective result at nmax = 800. b) Same as (a), but for helium.

size of the radial box R0 and the basis size. For the present set of calculations in hydrogen we havefound that a box size of R0 = 500 was sufficient to reduce the error to be within 0.5%. For helium,
a box size of R0 = 500 was used as well. The size of the numerical basis in the present calculation
is determined by the maximum angular momentum (lmax = (N + 1)/2, where N is the order ofthe coefficient calculated), and the maximum principal quantum number (nmax ) chosen. For a gridsize of R0 = 500, we varied the principal quantum number up to nmax = 1400 (corresponding toan energy of approximately 38.5 atomic units) for our perturbative calculations.
In Figure 4.4(a) we show results for χ(N ) as a function of nmax for R0 = 500 at ω = 0.056.
One can clearly see that the results for each of the coefficients converge to a finite value for large
nmax . The cross on each line indicates that the result is within 1% of the corresponding value at
nmax = 800. As one would expect, each higher order term requires additional states to be includedin the calculation in order to reach convergence. Furthermore, the dashed line at nmax = 19
indicates the separation of bound from continuum states in our basis set. For each of the nonlinear
coefficients (χ(N ) with N ≥ 3) we observe that the contributions from the bound states to the62coefficients are positive while those from the continuum states are negative. Thus, in order to reachconvergence obviously both bound and continuum states have to be included.

104
103
102
101
100
10−1
10−2
10−3
10−4
10−5
10−6
10−7
10−8
10−9
1012

Breakdown of the perturbative power series expansion104
χ(5) I
χ(3)102

χ(7) I 2
χ(3)χ(9) I 3χ(3)χ(11) I 4
χ(3)
P

1013

1014

intensity (W/cm2 )

(a)

1015

(N−1)

χ(N) I 2
χ(3) I

χ(N ) I (N −3)/2 /χ(3)χ(N ) I (N −3)/2 /χ(3)4.2.3.2

100
10−2
10−4
10−6
10−8
1012

1013

1014

1015

1016

1017

intensity (W/cm2 )(b)

Figure 4.5: a) The perturbative nonlinear susceptibility terms χ(N ) I (N −1)/2 scaled by χ(3) I as afunction of the intensity of a uniform electric field at 800 nm for the hydrogen atom. b) The samefor helium

In order to be applicable, the power series expansion, Equation (4.2), must converge. Since
we cannot determine all the terms in the expansion, we estimate a region of applicability by
investigating the relative contributions from the low order nonlinear terms in the expansion. To
this end, we have determined the ratios of the higher order nonlinear terms (N ≥ 5) to the lowest
order nonlinear term, i.e. χ(N ) I (N −1)/2 /(χ(3) I) as a function of the intensity I = |E|2 . The resultsfor hydrogen are shown in Figure 4.5(a) along with the sum of all higher-order terms calculatedwith respect to the third-order term (solid line). For low intensities (≈ 1012 W/cm2 ), each higherorder nonlinear term is much smaller than the lowest-order nonlinear term in the expansion. We
therefore can expect that the power series does converge in this intensity regime. On the other
hand, at high intensities (≈ 1014 − 1015 W/cm2 ) the trend is reversed, higher-order terms are larger
in magnitude than the lowest-order term, indicating that the series may diverge. For helium, this

XH = 0.01
XH = 0.1
XHe = 0.01
XHe = 0.1

N=5
2.6
26
26
260

N=7
0.91
9.1
8.0
80

N=9
0.63
6.3
7.2
73

N=11
0.29
2.9
3.2
32

63

Table 4.1: Intensities (in 1012 W/cm2 ) at which the ratio XH = χ(N ) I/χ(N −2) equals either 0.01
or 0.1 for hydrogen. XHe corresponds to the same for helium.

trend does not reverse until approximately 1015 − 1016 Wcm−2 , an order of magnitude higher thanfor hydrogen. This can be explained by the slow increase of the successive coefficients for hydrogenmentioned earlier.
In order to further narrow down the intensity regime at which the convergence of the power
series expansion becomes questionable, we note that the ratio of the sum of all higher order terms
calculated to the lowest-order nonlinear term exceeds 10% at about 3×1013 W/cm2 for the hydrogenatom. From the results in Table 4.1 we further see that in this intensity regime the ratio of successive
terms in the perturbative power series expansion does exceed 10%, which may be considered asanother indication for the inapplicability of the series expansion. This behavior is repeated for the
helium atom, but at a higher intensity. The ratio of the sum of all higher order terms calculated
to the lowest-order nonlinear terms exceeds 10% at about 3 × 1014 W/cm2 for the helium atom.
4.2.3.3

Comparison of results of perturbative and ab-initio calculations

Ab-initio calculations of the susceptibility are performed for laser pulses at a central wavelength of 800 nm having a sin2 -envelope and a total pulse length consisting of a finite numberof electric field cycles. In order to compare our results for the power series expansion, which are
determined for a uniform electric field, to those of the ab-initio calculations, we therefore need toaccount for the finite length and envelope of the pulse. To this end, we have averaged each term of
the power series expansion over the intensity distribution of the pulse.
The comparison of the results of the perturbative calculations and those of the ab-initio
calculations for a 10-cycle pulse and a 30-cycle pulse for hydrogen (Figure 4.6(a)) shows excellent

64
agreement at low intensities. In this regime of intensities all three results increase linearly withintensity, in agreement with the expectation that the lowest-order nonlinear term is dominantat these intensities. At higher intensities however the results of the two calculations reveal an
opposite trend. While the ab-initio results show, in good quantitative agreement with earlier
predictions [75, 76], and seen in Figure 4.1, a negative slope for the susceptibility as a functionof intensity, the sum of the first five terms in the power series expansion has a positive slope, as
expected from the positive sign of each term in the power series determined here. The 30 cycle
pulse shows a larger deviation from the perturbative calculations as the intensity increases, but
starts the deviation at the same general intensity as the 10 cycle pulse. Since the saturation of thesusceptibility is related to the ionization of the underlying media, we expect longer pulses to ionize
more, and thus deviate more strongly from the first order nonlinear susceptibility.
The peak intensity at which the discrepancy between perturbative and ab-initio results equals
1% is given in Table 4.2 for different pulse lengths at 800 nm. The result is slightly larger than our

5.5

0.90
0.85
0.80
ǫ0 χ(ω)/ρǫ0 χ(ω)/ρ

5.0

4.5

4.0

3.5
0.0

0.75
0.70
0.65

perturbative (800nm)
ab-initio (30 cycles, 800nm)
ab-initio (10 cycles, 800nm)

0.2

0.4

0.6

0.55

0.8

intensity (1014 W/cm2 )(a)perturbative
ab-initio

0.60

1.0

0

1

2

3

4

5

6

intensity (10 W/cm )
14

2

(b)

Figure 4.6: a) Comparison of results of perturbative (black line) and ab-initio calculations for thesusceptibility of hydrogen as a function of peak intensity of a 10-cycle and 30-cycle laser pulse at
800 nm. b) Comparison of results of perturbative and ab-initio calculations for the susceptibilityof helium as a function of peak intensity of a 10-cycle laser pulse at 800 nm.

pulse length
peak intensity

4
4.7

10
4.2

16
3.9

65

Table 4.2: Peak intensities (in 1013 W/cm2 ) at which the relative error between the results forthe ab-initio susceptiblity and the perturbative susceptibility equals 1% for laser pulses of different
pulse lengths (given in field cycles) at 800 nm for hydrogen.

estimates of the applicability limit of the perturbative series expansion in the previous subsection,
since the latter one has been obtained for a uniform field without ionization. However, we can stillconclude that the negative slope of the susceptibility is a result of the nonperturbative interaction
between the hydrogen atom and the intense laser pulse. Thus, an explanation of this featureusing higher-order expansion terms (e.g., HOKE) appears to be questionable and not applicable.
This also implies that a determination of higher order terms from ab-initio results in the intensity
regime above 2 − 4 × 1013 W/cm2 is meaningless. This is in agreement with the results of Köhler et
al. [29], who were not able to extract consistent HOKE terms from results of ab-initio simulationsin this intensity regime. Instead, if necessary, the atom-field dynamics has to be treated using anonperturbative theoretical description of the polarization and susceptibility.
For helium (Figure 4.6(b)), we see the same general trend: good agreement between the
results of the perturbative and ab-initio calculations at low intensities, and opposite trends athigh intensities, however for helium the higher order terms have a visibly smaller effect on theperturbative susceptibility. The third order susceptibility term dominates the perturbative solution
until past the divergence between the perturbative and ab-initio calculations.
The disagreement between the ab-initio and perturbative calculations for hydrogen and helium make clear the expected deviations for different ionization energies. We have further studied
this relationship between ionization energy and onset of the breakdown of the perturbative series
by using a Coulomb potential with different nucleus charges:
V (r) = −

Zeff
,
r

(4.11)

where we may vary Zeff , the nuclear charge, continuously. We show in Figure 4.7 the intensities
where a divergence occurred for 1% and 10% relative error as calculated in Table 4.2. We note that

66
8

2

intensity (W/cm )7
6

10% error
1% error

5
4
3
2
1
0
0.5

0.6 0.7 0.8 0.9 1.0 1.1
ground state energy (a.u.)

1.2

Figure 4.7: The intensity where the relative error between the perturbative and ab-initio calculations diverge by 1% and 10%, as a function of the ground state energy of atoms with nuclear charge
1, 1.25, and 1.5. The ab-initio calculations were done for a 10 cycle pulse with center frequency
corresponding to 800 nm. Convergence was checked in the same manner as for previous results.

the ionization energy correlates with the breakdown of the perturbation theory as expected, giventhat a higher ionization energy leads to smaller perturbative terms.

4.3

Low-order harmonic generation
The comparison between the perturbative series and the nonperturbative susceptibilities let us

identify an intensity regime for the breakdown of the perturbative series for hydrogen and helium.
We will now analyze whether or not our conclusions regarding the applicability of perturbation
theory can be confirmed by a similar study of low order harmonic generation 2 .
Kolesik et al. proposed a litmus test to determine if the higher order susceptibilities accurately
describe the physics of the propagation of high intensity light through a gas [33]. They did two
propagation calculations: one with higher order susceptibility terms, and one where they truncatedthe perturbative series after the first nonlinear susceptibility term. They then examined the thirdand fifth harmonic yields from both calculations. They found that the ratio of the fifth to third
2
Part of the results presented in this section are published in A. Spott et al., Transition from perturbative tononperturbative interaction in low-order-harmonic generation, Physical Review A 91, 023402 (2015).

67
harmonic yields was roughly 0.1 for the model that included the full perturbative series, while it
was significantly less than 10−3 for the model that truncated the series after the third harmonic
term, at the saturation intensity of the susceptibility (roughly 1014 Wcm−2 for the system they
were examining). They proposed that a measurement of this ratio might lead to a clear ruling on
the importance of higher order susceptibilities to the propagation of high intensity light.
There have been two experiments that have attempted to discern which model is correct.
The first, done by the same group that proposed the test, found that for femtosecond pulses withintensities measured in the 10s of TW/cm2 and a central wavelength of 2.2 µm the ratio of theyields between the third and fifth harmonics was 2 × 10−4 [35]. This was consistent with models
that did not have any higher order (greater than the third order) susceptibility terms, appearing
to confirm their hypothesis that the higher order Kerr terms did not play a role in high intensity,
short pulse propagation.
However, around the same time, another group performed similar experiments, with roughly
the same intensities and pulse duration, but a 1500 nm central wavelength [34]. This experiment
found that the ratio of the fifth over the third harmonic was closer to 0.1 at the saturation intensity,
giving credence to the higher order susceptibilities being necessary for properly describing the
propagation of high intensity laser pulses.
Using our method already developed, we believe that we can contribute to explaining thediscrepancy between these two experiments, and further expand our analysis from Section 4.2.3
towards longer wavelengths.

4.4

Perturbative and ab-initio calculations of low order harmonics
We calculated the susceptibilities for the low order harmonics using Equations (4.9) and(4.10), this time examining the 3ω and 5ω dipole moment contributions.
In Figure 4.8 we present an example for an ab-initio calculation of a low-order harmonic
spectrum generated at a driver wavelength of 1600 nm, a peak intensity of 5 × 1013 W/cm2 and apulse duration of 10 cycles. The results have been determined for a box size of R0 = 1000 a.u., a

68

10

−1

relative error

10−3
power (a.u.)

10−5
10−7
10

−9

103
101
10−1
10−3
10−5

0

5

10

15

20

10−11
10−13
10−15

0

5

10
harmonic order

15

20

Figure 4.8: Results of ab-initio numerical calculations for a low-order harmonic spectrum generated
by a driver laser pulse at a central wavelength of 1600 nm, a peak intensity of 5 × 1013 W/cm2 anda pulse length of 10 cycles. The inset shows the relative error between calculations using radial box
sizes of Rmax = 500 and Rmax = 1000.

time step of ∆t = 0.05 a.u. and a maximum principal quantum number nmax = 2000 as well as
lmax = 70. The convergence of the results with respect to the size of the radial box is shown by
the relative error between the results for box sizes of R0 = 500 and R0 = 1000. Please note thatthe minima in the error correspond to the peaks in the harmonics.
As for the ab-initio calculations, we performed test calculations to ensure that the resultsof our calculations for the terms in the perturbation expansion of the susceptibility for low-orderharmonic generation are converged with respect to the size of the box R0 and the size of the basis
set nmax . We note that the maximum angular momentum is determined by lmax = (N +1)/2, where
N is the order of the coefficient calculated. In general, we have found that a box size of R0 = 500 is(N )sufficient for the present purpose. In Figure 4.9 we show results for χnω for n = 1 (dashed-dottedlines), n = 3 (solid lines) and n = 5 (dotted lines) as a function of nmax for R0 = 500 at a laser
wavelength of 1600 nm. In general, we observe that the contributions from the bound states are
positive, reflected in the increase of the susceptibilities for low nmax , and those from the continuum

69

10

26

|ǫ0 χ(N ) (mω)/ρ| (a.u.)

1022
1018
1014
10

N = 11

10

N =9
N =7

106

N =5
N =3
N =1

102

10−2

0

200

400

600
800
nmax

1000 1200 1400

(N )Figure 4.9: Results for perturbative power series coefficients χnω for n = 1 (dashed-dotted lines),n = 3 (solid lines) and n = 5 (dotted lines) as a function of nmax .

states are negative, corresponding to the decrease of the susceptibilities for higher nmax .

4.4.1

Intensity dependence of low-order harmonic generation

We have applied both calculation methods to investigate the intensity dependence of low orderharmonics and the transition from perturbative to non-perturbative interaction in hydrogen. In
Figure 4.10 we present the integrated power of the (a) 1st, (b) 3rd and (c) 5th harmonic as a functionof peak laser intensity at a central wavelength of 1600 nm and a pulse length of 10 cycles, as obtained
from our ab-initio calculations. For these results we have calculated the harmonic spectrum andintegrated the signal for the power of the nth harmonic over the energy range [(n − 1)ω, (n + 1)ω].
We compare the results of our numerical calculations with the power law I n , which is expected for a
perturbative n-photon process. The predictions from the power law were matched to the numericalresults at low intensities. The inset in each of the panels shows the relative error between the
ab-initio results and the power law predictions with respect to the ab-initio results.
The results show that in the intensity regime between 1012 Wcm−2 and a few times 1013 Wcm−2
the ab-initio results start to deviate from the respective power law. This is an indication of the

70

10

0.00
−0.05
−0.10

−5

1012

1013

10−6
(a)
integrated power (a.u.)

1012
10−7
10−8
10−9
10−10
10−11
10−12
10−13
10−14
−8

1012

10
10−10
10−12
10−14
10−16 11
10
10−18
10−20
10−22
10−24 11
10

1013
0.00
−0.05
−0.10

1013

(b)
1013

1012

1012

1013

0.00
−0.05
−0.10
(c)

1012

1013

intensity (W/cm2 )Figure 4.10: Results of ab-initio calculations for the integrated harmonic power (solid circles withsolid lines) for the (a) 1st, (b) 3rd and (c) 5th harmonic as a function of the peak laser intensityof a laser pulse of 10 cycles at a wavelength of 1600 nm. The numerical results are compared toa perturbative I n power law fit, which is matched to the ab-initio results at the lowest intensity.
The insets show the relative error between ab-initio results and power law predictions with respect
to the ab-initio results.

transition from a perturbative to a non-perturbative electron-field interaction. These results are inagreement with previous results that the onset of other non-perturbative phenomena, e.g., above
threshold ionization [3] and high-order harmonic generation [77, 78], occurs in the same intensity
regime.

71

Figure 4.11: Results for the ratio of higher-order terms to the lowest-order nonlinear term in theperturbative series expansion for (a) χω , (b) χ3ω and (c) χ5ω . Also shown is the ratio of the sumof all higher order terms calculated with respect to the lowest-order term (solid lines).

Based on the ab-initio results, we expect that the perturbative power series expansion of thesusceptibility corresponding to the process of low-order harmonic generation should break down inthis intensity range as well. In order to test this expectation, we have calculated the first few terms

72
in the expansion for χω , χ3ω and χ5ω at 1600 nm. We study the relative contribution of higherorder terms in the expansion by presenting their ratios with respect to the lowest-order nonlinear
term in Figure 4.11. Also, shown is the ratio of the sum of all higher-order terms calculated to
the lowest-order term. From the results we observe the same behavior, as previously reported forthe electrical susceptibility χω at shorter wavelength [79], namely each higher order term is muchsmaller than the lowest order term, indicating the convergence of the corresponding power series
expansion, at the lowest intensities studied. On the other hand, the breakdown of the series at
the highest intensities is obvious as well, since the contributions of the higher order terms exceedthat of the lowest order term. We further note that in each case the sum of the calculated higherorder terms reaches about 10% of the lowest-order term for intensities in the range of 1 × 1013 to2 × 1013 W/cm2 . Therefore, this limit can be considered as an indication for the breakdown of a
perturbative series expansion in strong-field processes for the hydrogen atom [23].

4.5Summary
We have performed calculations of the susceptibility of atomic hydrogen and atomic helium

using a numerical basis state method, via a perturbative power series expansion as well as theab-initio results for the field induced polarization. The results of our perturbative calculations
indicate a breakdown of the series expansion at intensities in the lower range of 1013 Wcm−2 forhydrogen and 1014 Wcm−2 for helium. Comparison with the ab-initio results shows a discrepancy
in the trend in the same intensity regime. We therefore conclude that the negative slope of thesusceptibility of atomic hydrogen and helium at higher intensities cannot be explained by the higherorder terms of a perturbative power series expansion. Instead, the feature has to be considered asa signature of the nonperturbative interaction between the atom and the field.
Our results for low-order harmonic generation from both ab-initio as well as perturbative calculations show the same onset of a transition from perturbative to non-perturbative interaction between the atom and the field as the results for the electrical susceptibility χω at shorter wavelength.
It is therefore not surprising that previous studies on low-order harmonic yields [26, 34, 35, 80, 81]

73
did not help in resolving the question about the significance of higher-order Kerr effects in thefilamentation of short higher-power laser pulses in gaseous media. In contrast, we conclude thatany deviations from the predictions of the lowest-order perturbation theory for the polarization(and other observables) should be interpreted as a signature for the non-perturbative character ofthe electron-field interaction. In particular, our results also show that a quantitative analysis ofstrong-field below-threshold harmonic generation requires a nonperturbative theoretical approach,as e.g. introduced in [82].
While we are not able to prove a general statement about the breakdown of the perturbative
susceptibility, combined with the previous results showing ATI peaks in xeon [3], and high harmonic
generation [17] at similar wavelengths, we believe that a strong argument has been made for the
inapplicability of perturbative descriptions of strong field processes around these intensities, andfor these wavelengths.

Chapter 5Time Dependent Susceptibility

In this Chapter, we introduce a method to calculate the time dependent susceptibility using
the short time Fourier transform. First, we analyze various windowing functions for determining
the short time Fourier transform and settle on the Dolph-Chebyshev family of windowing functions
for our time dependent susceptibility calculations.
We then examine the time dependent susceptibility of helium at intensities near the breakdown of the perturbative susceptibility found in the previous Chapter. Using this method, we
examine several features that show up at these intensities and explain them via the dynamics ofthe field-free bound state and continuum state contributions to the susceptibility. From there,
we further study the emergence of the non-perturbative multi-photon ionization evident in the
population growth and the susceptibility of helium as we increase the intensity.

5.1

Introduction
High harmonic generation (HHG) is a process for generating extreme ultraviolet (10 ∼ 100eV), soft x-ray (100 ∼ 1000 eV) and even hard x-ray (> 1 keV) coherent photons from much
longer wavelength driving fields. The success in generating bright beams is largely a combinationof microscopic and macroscopic phenomena. On the microscopic level, HHG is a process in which
an electron in a single atom is driven by the external field and generates radiation of large multiplesof the frequency of the driving field, as discussed in section 1.2 previously. On the macroscopic
level, a coherent constructive interference of these high harmonic emissions from a large number of75
atoms needs to be achieved in order to generate a bright laser beam [83–85]. Ensuring that thisinterference is constructive instead of destructive in the direction of the propagation of the driving
field is mainly done by matching the phase velocity of the driving light field with the phase velocityof the generated light. We note that geometric effect as well as other effects, such as group velocity
matching [86, 87] can play a role as well, depending on the experimental setup and the wavelengthof the driving field.
In general, this phase matching relationship can be written as (e.g. [88]):∆k(t) = qkω (t) − kqω (t)(5.1)

where q is the harmonic number, kω (t), and kqω (t) are the magnitudes of the wavevector of thedriving field and the qth harmonic, respectively, and ∆k is the phase mismatch between them. In
terms of the index of refraction, this can be written as∆p(t) = n(ω; t) − n(qω; t) + ∆g(t)(5.2)

where ∆p(t) is the phase velocity mismatch, n(ω; t) is the index of refraction of the driving field,n(qω; t) is the index of refraction of the qth harmonic, and ∆g(t) is a geometry induced factor forthe beam parameters. In this form, there is an explicit relation between the time dependent index
of refraction and the phase velocity mismatch.
Typically, the phase matching is determined by the n(ω; t) term, since the index of refractionof the higher harmonics (n(qω; t), for large q), is well approximated by 1. Thus, the phase velocity
mismatch is described in terms of the index of refraction of the fundamental field, which can be
split up into two contributions (e.g. [87]):n(ω; t) ≈ natom (ω; t) + nion (ω; t)(5.3)

The term natom (ω; t) is the contribution from the neutral atoms or molecules. This term is typically
positive, and will approach zero as the neutral atoms are ionized during the pulse, and thus thenumber density of neutral atoms decreases. The second term nion (ω; t) is the index of refraction
due to the ionized electrons. This term is negative, and will increase in magnitude as more atoms

76
are getting ionized. HHG is often observed in low density gas and the nonlinear terms in thesecalculations are typically ignored, leading to the index of refraction being approximated by:(1)ni (ω, t) ≈ 1 + 21 χi (ω)ηi (t)

(5.4)

(1)

where i represents either the atoms or the ionized electrons, χi (ω) is the linear susceptibility, andηi (t) is the time dependent density of the atoms or ionized electrons.
Recently, however, two results have challenged this conventional approximation of a linear susceptibility: the generation of extremely high harmonics from very powerful mid-infraredlasers [89], and the generation of bright harmonics in multi-ionized noble gases using powerful UV
lasers [90]. In the first case, the generation of x-rays from mid-infrared lasers is improved by filamentation using very high pressure gas. From Chapter 1, we know that filamentation occurs in theregime where nonlinear effects can no longer be neglected and even not described perturbatively.
We may therefore expect that a more powerful ab-initio non-perturbative method of determining
the time dependent index of refraction might be useful. In the second case, the generation of brightharmonics by powerful UV laser sources requires intensities where multiple ionization occurs, which
is likely beyond the saturation point of the nonlinear refractive index as a function of intensity. As
we have previously shown in Chapter 4, this intensity regime marks a breakdown in the accuracyof the conventionally used perturbative series for the refractive index.
Examining the nonlinear response near the ionization threshold as a function of time is also
another useful test towards further understanding the process of filamentation and the limits ofperturbation theory. In investigating this question, Wahlstrand et al. [91–93] developed a techniquefor observing the time dependent response of a thin gas to a strong pump pulse by measuring
the phase delay of a weak supercontinuum probe pulse coincident with the pump pulse. The
supercontinuum probe pulse is chirped, in order to read out the phase shift as a function of timeand determine a time dependent response.
The experimental results performed by Wahlstrand et al. for neon below and above the onsetof ionization are presented in Figure 5.1. It shows both the transverse extent of the pump pulse,

PRL 109, 113904 (2012)PHYSICAL REVIEW LETTERS
77

obtained from the sernegative in the inten
below the threshol&80 TW=cm2 . At ve
reasonable agreemen
experimental results
Equation (1) for
absorption spectrum
and !0 , which availa
optical fields cannextended an ‘‘avera
by Sheik-Bahae etto arbitrary ord""ð!1 ; !2 ; . . . ; !K Þ
sorption of K pump p
photon of frequency
written as

1 (color
online).
Probephase shift"#ðx;measuredin
Figure 5.1: TheFIG.
phase
shift ∆Φin radians
as a functionof timeand tÞ
transverse
positionin a gas ofneon
usingaresingle-shotsupercontinuum
spectralneon. On the leftshown
results observedat an intensity
below interferometry,
the ionization threshold, while
for right
an intensity
below (left)
above (right)
the Note
onsettheofnegative phase
the results on theare at an intensity
above andthe ionizationthreshold.
ionization.
shift from the ionizedelectrons present on the right, after the probe pulse has ended (From[91]).!
where

""ð!;

avg;K ¼ ðK!

tation with a single aapproach, but examin
been reviewed by Hutchings, Sheik-Bahae et al. [10], andIn any considerati
has been used to calculate the real nonlinear index ofand the phase shift as a function of time. The phase shift is directly related to the shift in thesorption mechanism
refraction (n2 ) of semiconductors from two-photon absorpaware
index of refraction
due to transitions
the interactionof the pumppulseeffectwith thegas.InHowever,resultsof the transiti
tion,∆nRaman
andac Stark[11].
this thenonperturbative (tunnanalysis, they related the real and imaginary parts of themust be carefully interpreted, as each time is also associated with a corresponding frequency in theof this transition is t
response quadratic in the field strength.
supercontinuum signal,
resultingin ∆n(ω(t),with ω andt coupled. absorption$ ¼ ð2%p Þ1=2 !=E0 , w
To apply
the general
KKt),relation
to nonlinear
and the nonlinear refractive index [10], one considers thea gas atom, E0 is the
In order to explain the experimental results, we use the theoretical ab-initio approach devellinear response to a weak probe field at frequency ! of theused. While one is sa
oped in the previous
Chapters
to determine
a time dependent
susceptibility
analyze theit. Intunnelingthesystemconsistingof the matterand a strong
optical andpump
limit fo0
field Eð! Þ. The real refractive index change "n for aargued that $ < 1 isrest of this Chapter, we will first give an overview of the method of calculating the time dependent
probe field at frequency ! is thenand $ > 1 to be in thsusceptibility, and then discuss some results for the helium atom at0 different laser parameters.
when considering thc Z 1 ""ð!; Eð! ÞÞ0erably more careful.
"nð!; Eð! ÞÞ ¼ P
d!;
(1)
!
!2 $ !20
In their calculation
k-photon ionization
where P denotes principal value. Here, ""ð!; Eð!0 ÞÞ is
MPI limit of a recentl
the change in absorption coefficient of a probe beam aton the theory of Pe
frequency ! in the system composed of the matter and aexpression R (Eq. 6 opump field Eð!0 Þ. Below we are interested only in the$. It handles the MPI
index change for nearly degenerate fields, so we set
accurate interpolatio!0 ¼ ! after performing the integration [10]. For thea single active elec
Kerr and ionization processes discussed in this Letter,
In Fig. 2(a), we plot
where the response time is essentially instantaneous and

78

5.2

Theoretical Method
A non-perturbative analysis of the time dependent susceptibility requires the calculation of

the time dependent dipole moment beyond the conventional approach of perturbation theory. Thus
we use our ab-initio method described in Chapter 3 to determine the dipole moment. As a brief
reminder, we find the time dependent dipole moment as the expectation value of the dipole operatorµ̂ with respect to the time dependent wavefunction |Ψ(t)i as:µ(t) = hΨ(t)| µ̂ |Ψ(t)i(5.5)

where the dipole operator is defined as:µ̂ = −

Xn,l,n0 ,l0|ψn,l i hψn,l |z|ψn0 ,l0 i hψn0 ,l0 | ,(5.6)

and |ψn,l i = Rn,l (r)Yl (r̂) are the eigenstates of the field-free Hamiltonian. The wavefunction |Ψ(t)iis the solution to the time dependent Schrödinger equation with Hamiltonian:H = H0 + µ̂E(t)H0 =∇2+ V (r).
2(5.7)

Here, we assume that V (r) is a spherically symmetric, single active electron potential, and H0 is
the field-free Hamiltonian.

5.2.1

Time dependent susceptibility

The time independent susceptibility is given by:χ(ω) =

F {µ(t)} (ω)
F {E(t)} (ω)(5.8)

where F {x(t)} is the Fourier transform of x(t). If one multiplies the dipole moment µ(t) and theelectric field E(t) by a windowing function before taking the Fourier Transform, one obtains a timedependent susceptibility:χ(ω; τ ) =F {µ(t)w(t − τ )} (ω),
F {E(t)w(t − τ )} (ω)

(5.9)

79
where the time τ (corresponding to the center of the windowing function) is now a parameter tothe susceptibility, and w(t) is a windowing function with limited support, i.e. it is only non-zeroover some finite range [−T /2, T /2].
The windowing of the Fourier transform is known as the short time Fourier transform [94]. By
limiting the non-zero part of the time domain function one defines a time associated with the Fouriertransform. By moving the center of the windowing function, one gets a Fourier transform for eachpoint in time. When one applies this concept to the susceptibility, one obtains the time dependent
susceptibility given in Equation (5.9). As the properties of the short time Fourier transform depend
strongly on the windowing function that is being used, we will first discuss some of the propertiesof the windowing function, before presenting the results for the time dependent susceptibility.

5.2.2

Windowing functions

The properties of the windowing function are important to get an accurate time dependent
susceptibility. A windowing function that is too wide in frequency will not allow us to separate
out the driving frequency behavior from the harmonics (primarily the third and fifth), or the low
frequency drift of ionized wavepackets. On the other hand, one that is too wide in time will notgive the necessary time resolution and will effectively smooth out the short time behavior.
A windowing function w(t − t0 ) is any function that is non-zero only for some limited width intime, and is symmetric about some center time t−t0 = 0. The simplest windowing function is known
as the boxcar or rectangular window [95]. This function is 1 everywhere within the window size.
Applying it to a function f (t) is equivalent to just taking the Fourier transform of the overlappedw(t − t0 )f (t).
The boxcar window, however, is a very poor window for many applications. The sharp cutoff
at the edges of the boxcar window creates large sidelobes when taking the Fourier transform. The
Fourier transform of such a square pulse is the sinc function (= sin(ω)/ω), which results in a 1/ω
falloff of the sidelobes. This can be seen in Figure 5.2, which shows the amplitude of the discrete
Fourier transform on a log scale. One can clearly see the sidelobes besides the central peak and

80
1.04

100

100

10 1

10 1

power

1.00

power

amplitude

1.02

0.98
0.96
0 5 10 15 20 25 30 35 40 45 50
samples

20

10

0
10
freq. bins

20

10 8 6 4 2 0 2 4 6 8 10
freq. bins

Figure 5.2: The boxcar windowing function (a) in the time domain, (b) the frequency domain (logscale), and (c) a close up view of the frequency domain (log scale) with the harmonic orders denoted
for a single cycle window. Note that the attenuation for the first sidelobe is very poor (less than
10db)

their considerable height. The x-axis is presented in units of the discrete Fourier transform (DFT)bins – the smallest frequency difference that is representable by a DFT of the same size as thewindow.
Since multiplying a function f (t) by a windowing function in the time domain is equivalent
to taking the convolution of the windowing function and fˆ(ω) in the Fourier domain, the above example shows that the frequency response of the windowing function is important for understanding
the behavior of the window function on the signal. There are a couple of criteria that are importantwhen choosing a windowing function. For our purpose, we wish to optimize with respect to thefollowing features:(1) Small influence of the contributions from odd harmonics of the driving frequency. The 3rd(and 5th) harmonic should not affect the calculated susceptibility at the driving frequency,
thus the sidelobes at these frequencies should be suppressed.
(2) Small interference between the contributions from the driving frequency and the very
low frequency part of the dipole moment. During strong-field interactions, an electron
wavepacket that is ionized will begin to drift with a velocity that is proportional to thevector potential at the time of ionization. This will show up as a component at a very low

81frequency to the susceptibility which should not have an effect on the calculated results atthe driving frequency in the calculations.
(3) Small window size in time, in order to examine the processes with good time resolution.
The sidelobe attenuation defines the effect that the higher order harmonics and the DC componentof the field will have on the signal. If the sidelobes of the windowing function are not suppressedenough, then there will be some contribution overlapping with the time dependent susceptibility at
the frequency of interest. The central peak bandwidth is also important for controlling how muchof the signal from other frequencies effects the result at the central frequency. If the bandwidth is
too large, then the calculated signal at the primary frequency will include contributions from theDC or even third harmonic components. These properties must be taken into account in order to
optimize the result with respect to the first and second points above. The central bandwidth of any
windowing function – as well as the sidelobe suppression in many cases – is inversely proportional
to the size of the windowing function in time. This means that one must balance a narrow central
peak with high sidelobe suppression in order to avoid the influence of unwanted signals from otherfrequency components, with a large central bandwidth in order to have better time resolution.
There are a few windowing functions that have the required behavior, which will be discussed
below. The central peak of a window in the “flattop” family of windowing functions is very flat,
giving it its name. Windows in this family are typically cosine windows (windows defined as a sumof scaled cosines), and are optimized such that the derivative in the frequency domain at ω = 0 is
zero. This property is valuable when the sampling time and the periodicity of the signal are notwell aligned since the frequency of interest will be captured even though it is not well matched withthe frequency bins.
A flattop windowing function from [95] that has been simultaneously optimized for large
sidelobe attenuation and a flat top is shown in the time and frequency domains in Figure 5.3. Thiswindow is called a 5th order window, consisting of five cosine terms, and as such the first sidelobe

82
1.0

power

0.6
0.4

100

10 1

10 1

10 2

10 2

10 3

10 3

10 4

10 4

power

amplitude

0.8

100

0.2
0.0
0 5 10 15 20 25 30 35 40 45 50
samples

20

10

0
10
freq. bins

20

10 8

6

4

2 0 2
freq. bins

4

6

8 10

Figure 5.3: Flattop window given in Equation 5.10, in the time domain (a), the frequency domain(b), and with a closeup of the frequency domain (c). Note that the wide central peak covers 5
frequency bins.

will be the 5th frequency bin from the center. This flattop window is given approximately by [95]:w(n) = 0.216−0.417 cos










2πn
4πn
6πn
8πn
+0.277 cos−0.084 cos+0.007 cos
,
N −1
N −1
N −1
N −1
(5.10)

where N is the total size of the window in samples and n is the number of the current point insamples.
A window which has been optimized for the highest attenuation in the first sidelobe (thesidelobe nearest to the center peak) provides another interesting option. The Hamming window [96]is a second order window with such a property,
wHamming (n) = 0.54 − 0.46 cos



2πn
N −1



,

(5.11)

while the Blackman-Harris window [96] is a third order window of this kind
wBlackman-Harris (n) = 0.359−0.488 cos








2πn
4πn
6πn
+0.141 cos−0.012 cos. (5.12)
N −1
N −1
N −1

The first sidelobe of the Hamming window has an attenuation of roughly a factor of 43db, while
the Blackhman-Harris window has an attenuation of roughly a factor of 65db.
Finally, the Dolph-Chebyshev window [97], shown in Figure 5.4, is a window that is optimized
to have the narrowest central bandwidth for a given constant attenuation in the sidelobes. The

83
1.0

power

0.6
0.4
0.2
0.0

0 5 10 15 20 25 30 35 40 45 50
samples

100

10 1

10 1

10 2

10 2

10 3

10 3

10 4

10 4
20

10

0
10
freq. bins

20

10 8

6

4

2 0 2
freq. bins

4

6

power

amplitude

0.8

100

8 10

Figure 5.4: The Dolph-Chebyshev window with an attenuation set to 80db in the (a) time, (b)frequency and (c) a closeup of the frequency domain

Dolph-Chebyshev window is defined in the Fourier domain, rather than the time domain:]cos N cos−1 [β cos πkN
F {w(n)} (k) =−1cosh N cosh (β)
 A 
1−1
10 20 ,
β = coshcoshN(5.13)(5.14)

where A is the attenuation in decibels, and k is the frequency component.
None of these windows have a central lobe that has a full width of less than four frequency
bins. This necessitates our window encompassing more than one period of our target frequency,
in order to avoid leaking from the DC component. In the rest of this Chapter we use the DolphChebyshev window with a width of three periods of the frequency of interest. Our test calculations
have shown that it is a good compromise between attenuation of components that we are notinterested in, and a reasonably good time resolution.

5.2.3

Segmentation of time dependent susceptibility

Our use of a basis state method makes any observables accessible for separation in contributions from the eigenstates of the field-free Hamiltonian. Thus, we can split the time dependent
susceptibility into contributions arising from bound and continuum states in order to better understand the influence of these states on the total susceptibility. As pointed out before, in our method

84
the wavefunction in the energy basis is a coherent superposition of the field-free energy states:|Ψ(t)i =Xn,lcn,l,E (t) |φn,l,E i ,(5.15)

where the eigenstates |φn,l,E i are labeled by energy E as well as the principal quantum number n
and the angular quantum number l. Thus, by selectively setting cn,l,E to zero for some choice ofn, l and E, we can select the wavefunction to find a partial dipole moment.
We examine the wavefunction in two principle regions of the field-free spectrum: the boundstates consisting of all states with energy less than zero, and the continuum states consisting of allstates with energy greater than zero. This split then gives us three terms in the dipole moment:bound-boundcontinuum-continuum

z
}|
{ z
}|
{
†
hΨ(t)| µ̂ |Ψ(t)i = hΨ(t)| Pb µ̂Pb |Ψ(t)i + hΨ(t)| Pc† µ̂Pc |Ψ(t)i+zhbound-continuum}|hψ(t)| Pb† µ̂Pc |ψ(t)i + c.c.

i{

(5.16)

where Pb is the projection operator onto the bound states which sets cn,l,E = 0 for E > 0, and Pc is
the projection operator onto the continuum states. These three terms – referred to as bound-bound,bound-continuum and continuum-continuum contributions to the susceptibility – will be discussed
in the remainder of the Chapter.

5.3

Results
We have used the Short time Fourier transform to calculate the time dependent susceptibilityof Helium at five different intensities: 100, 150, 200, 250, and 300 TWcm−2 . These calculations
were performed for a 10 cycle full width at half max pulse at 800 nm. The results for the middle
16 cycles are shown in Figure 5.5, with a vertical dot-dashed line denoting the center of the electricfield of the pulse. We denoted the central cycle of the field by “0.” For the low intensities (100,150 and 200 TWcm−2 ), the susceptibilities adiabatically follow the envelope of the field intensity.
Above 200 TWcm−2 , however, the results deviates from this adiabatic behavior. At 250 TWcm−2
the peak of the susceptibility shifts towards the beginning of the pulse, and at 300 TWcm−2 it has

0.820

0.815

0.810

0.810

0.805

0.805

0.800

0.800

0.795

2

0.790
0.785

0.795

2

0.790

b)
−8

85

0.820

1

−6

−4

−2

0
cycles

2

4

6

8

0.0

2.5
5.0
cyclessusceptibility (a.u.)susceptibility (a.u.)

0.815

1

a)

7.5

0.785

Figure 5.5: (a) The time dependent susceptibility of helium for Gaussian pulses with a wavelength of
800 nm and a full-width at half-maximum of 10 cycles at different intensities (100 TWcm−2 (purple),
150 TWcm−2 (red), 200 TWcm−2 (green), 250 TWcm−2 (orange), and 300 TWcm−2 (blue)). The
center of the driving electric field pulse is denoted with a vertical dot-dashed line, and time isdenoted in cycles of the driving field. The arrows point out the shifted peak in the susceptibility(1) and the “kink” (2) at highest intensities. Note that both the 250 TWcm−2 pulse and the
300 TWcm−2 pulse peak before the driving field. (b) A detailed plot for 300 TWcm−2 .

moved by more than a full cycle of the driving field. This behavior is in qualitative agreementwith results presented in the supplemental materials of [91]. However, we note that it is evenmore obvious in the experimental results for krypton and xenon than for helium. After the peakin the susceptibility, the slope in the trailing edge for the 300 TWcm−2 pulse also deviates fromadiabatic behavior, exhibiting two regions over which the susceptibility follows an approximately
linear dependence on time: from -1 cycles to +5 cycles, and from +5 cycles on. This creates a“kink” at +5 cycles where the slope of the susceptibility changes from the strongly negative slope
between -1 and +5 cycles to the less negative slope after +5 cycles. In the rest of this section we
focus on the origin of these high intensity features – the “kink” and the shift of the peak of thepulse – within the context of our energy basis calculations. In particular, we analyze the field freebound and continuum state contributions to the susceptibility as a function of time.

5.3.1

Bound-bound contributions

The bound-bound contribution to the time dependent susceptibility (see Equation (5.16))
is shown in Figure 5.6(b), accompanied by the excited state population taken at the zeros of the

86
a)

1.2

b)

1.0
0.8

0.44

0.6
0.43

0.4

0.42
0.41

population (x10−5)

susceptibility (a.u.)

0.45

0.2
0.0
−8

−6

−4

−2

0
cycles

2

4

6

8

−8

−6

−4

−2

0
cycles24

68

Figure 5.6: a) Bound-bound contribution to the time dependent susceptibility, (hPb† µ̂Pb i (t)). It
follows the intensity envelope of the field for all intensities over the majority of the pulse. We notea small deviation at the end of the pulse for 300 TWcm−2 where the residual population in the
excited states results in a higher value of the susceptibility at the end of the pulse as comparedto the beginning of the pulse. (b) Population of the excited states taken at the zeros of thedriving field. Note the non-adiabatic contribution to the population for driving fields stronger than
150 TWcm−2 , along with the residual population at the end of the pulse for 300 TWcm−2 . For bothpanels: 100 TWcm−2 (purple), 150 TWcm−2 (red), 200 TWcm−2 (green), 250 TWcm−2 (orange),and 300 TWcm−2 (blue) pulses.

driving electric field (Figure 5.6(a)). The excited states are here defined as field free states withenergy less than zero, excluding the ground state. For the lowest intensity (100 TWcm−2 purple
curve), the population and the susceptibility adiabatically follow the envelope of the intensity, as
expected. For higher intensities, the population does not follow the electric field intensity smoothly.
Instead, the population consists of a dominant contribution driven by the perturbative interaction,
where the first order and third order terms of the susceptibility dominate, and a much smaller nonperturbative part, due to the population in excited states with larger l which are populated in the
trailing part of the pulse via transitions back from the continuum. We will further analyze this nonperturbative behavior in subsection 5.3.4. For the highest intensity considered (300 TWcm−2 ), theresidual population in the high l states results in a bound-bound contribution to the susceptibility
at the end of the pulse that is slightly larger than at the beginning of the pulse.
The susceptibility from the bound states is the largest contribution to the total susceptibility.

87
5
4

−0.2
−0.4

3

−0.6

2

−0.8
−1.0

1

−1.2

0
−8

−6

−4

−2

0
cycles

2

4

6

8

−8

−6

−4

−2

0
cycles2

4

6

population (x10−5)

susceptibility (x10−2 a.u.)

0.0

8

Figure 5.7: (a) Continuum-continuum contribution to the time dependent susceptibility(hPc† µ̂Pc i (t)). This contribution to the susceptibility is negative at high intensities (200 TWcm−2
and above). (b) The population of the continuum states taken at the zeroes of the driving field.
Note the exponential growth for 300 TWcm−2 , which starts slightly before the center of the pulse
and peaks slightly after. This is followed by a decline to a residual population at the end of the pulse.
For both panels: 100 TWcm−2 (purple), 150 TWcm−2 (red), 200 TWcm−2 (green), 250 TWcm−2(orange), and 300 TWcm−2 (blue) pulses.

It is therefore the reason that the overall shape of the susceptibility follows (almost adiabatically)
the envelope of the intensity. Indeed at low intensities, this contribution to the susceptibility
almost entirely defines the shape of the total susceptibility. For the deviations from field-following
behavior, specifically the features we pointed out earlier at higher intensities, we will now analyze
the other contributions in Equation (5.16).

5.3.2

Continuum state contributions

The continuum-continuum contribution to the time dependent susceptibility is shown in
Figure 5.7(a), accompanied by the population in the continuum states (with energy greater thanzero) taken at the zeros of the driving electric field (Figure 5.7(b)). At the lowest intensities, theshape of the susceptibility as a function of time is again adiabatic, and positive, though small.
At these intensities, there is no residual population in the continuum at the end of the pulse andthe population in the continuum is mainly transient and likely via perturbative interactions. This

88
agrees with our observation that the population in the continuum states exponentially decreases
with increasing energy. We show this in Figure 5.8, where the population is shown as a functionof the energy at the zero of the field at cycle +1. We also note that the population in these states
exhibits the same qualitative behavior as that in the excited states just below the threshold; namely
an exponential decrease in the population with increasing energy, and a positive contribution to
the susceptibility, as expected from perturbation theory [79]In contrast, at the highest intensities (250 and 300 TWcm−2 ), the continuum is strongly (nonperturbatively) populated between -2 and +1 cycles before relaxing back to a residual continuum
population at the end of the pulse. The exponential growth in population over a few cycles is a clear
indication of the deviation from the perturbative regime and correspondingly leads to similar characteristic changes in the susceptibility as a function of time. Figure 5.8 shows that the populationis almost equally distributed among the low energy continuum states, another strong indicationof a deviation from a perturbative interaction between the field and helium atom. We note that
this result agrees with our previous results in Chapter 4, where we have observed the transition

10−6
10−7

population

10−8
10−9
10−10
10−11
10−12
10−13
0

1

2

3
energy (a.u.)

4

5

6

Figure 5.8: The population in the continuum states as a function of energy at cycle 1 for
100 TWcm−2 (purple), 150 TWcm−2 (red), 200 TWcm−2 (green), 250 TWcm−2 (orange), and
300 TWcm−2 (blue) pulses. Note the exponential decrease in population for the 100 TWcm−2
pulse, while at the higher intensity pulses the distribution extends to higher energies.

89
from the perturbative to the non-perturbative regime at similar intensities. The population in thecontinuum from the multiphoton transitions to the continuum causes a negative contribution to
the susceptibility, in contrast to the signatures seen at low intensities. We also note that the dip in
the continuum-continuum contribution to the susceptibility at +7 cycles (Figure 5.7(a)) is matchedby the peak in the continuum population at +7 cycles in (Figure 5.7(b)).
These features in the contribution from the continuum states at the highest intensities can
explain the signatures we pointed out in the total susceptibility (Figure 5.5). The strong growthin the continuum population (and corresponding negative contribution from the continuum contributions to the susceptibility) shows up in the total susceptibility as a shift in the peak of thesusceptibility towards earlier times in the pulse. The peak in the population at around +1 cycles(and the corresponding minimum in the susceptibility) results in the linear regime between -1 and+5 cycles until the “kink” in the total susceptibility when the trapped population in the boundstates contributes significantly enough to the total susceptibility to cancel out part of the negative
contribution from the continuum states at approximately +5 cycles.

5.3.3

Cross term contributions

The contributions to the susceptibility from the cross terms (hΨ(t)| Pb† µ̂Pc |Ψ(t)i + c.c.) are
shown in Figure 5.9. Overall, the variation in this contribution is small and, hence, it is notimportant for our understanding of the dynamic features in the total susceptibility. We observe a
low field-contribution (at the beginning and end of the pulse), a field envelope following negative
contribution, and small deviations at higher intensities. While the overall offset is due to the firstorder linear term to the susceptibility, the field envelope following negative contributions result
from the third order term. The negative shape is likely due to a phase difference between the thirdorder perturbative population in the continuum and the ground state with which it is interacting(i.e. the hψ (3) (t)| µ̂ |ψ (0) (t)i term in the dipole moment). The deviations of this shape at higher
intensities are likely due to non-perturbative effects.

90
0.3745susceptibility (a.u.)

0.37400.37350.37300.37250.37200.3715
−8

−6

−4

−2

0
cycles24

68

Figure 5.9: Cross-term contributions to the susceptibility: (100 TWcm−2 (purple), 150 TWcm−2(red), 200 TWcm−2 (green), 250 TWcm−2 (orange), and 300 TWcm−2 (blue) pulses.) Note thelarge offset but subsequent small dip in the susceptibility. The dynamic variation here is much lessfor high intensities than for the bound states or the continuum states.

5.3.4

Populations and susceptibilities for high and low angular momentum

As mentioned previously, we attribute part of the (nonadiabatic) signatures in the totalsusceptibility at high intensities as due to the transition to a nonperturbative interaction of the fieldwith the target atom. In particular, the population of the continuum states clearly shows featuresof the nonperturbative interaction, e.g. the exponential growth in continuum state population overa short interaction time.
While it is difficult to separate perturbative and nonperturbative effects in our ab-initio
calculations, we noted another interesting signature of the transition between the corresponding
regimes. This is the population of excited bound states with angular momentum quantum number
larger than 1. To show this characteristic feature we present in Figure 5.10 the population in the
excited bound states, recorded at the zeroes of the field, as a function of interaction: (a) for allstates with l = 1 and (b) for all states with l > 1. At all intensities considered, we note that the
population in the p states (Figure 5.10(a)) mainly adiabatically follows the field envelope, while
the population in the higher l states appears in the trailing part of the pulse. Moreover, the latter
population is predominantly observed at higher intensities. This effect reminds of various previous

91
a)

b)

1.75

1.0

1.50

0.8

1.25
1.00

0.6

0.75

0.4

0.50

0.2

0.25

0.0

0.00
−8

−6

−4

−2

0
cycles

2

4

6

8

−8

−6

−4

−2

0
cycles2

4

6

population (x10−6)

population (x10−5)

1.2

8

Figure 5.10: Excited state population, for (a) l ≤ 1 and (b) l > 1. The lines represent resultsof calculations for 100 TWcm−2 (purple), 150 TWcm−2 (red), 200 TWcm−2 (green), 250 TWcm−2(orange), and 300 TWcm−2 (blue) pulses.

studies [69, 98, 99] in which the trapping of population in excited states via transition from thecontinuum in the trailing part of the pulse has been observed and analyzed. Here, our results showthat this excited population has significant contributions from states with l ≥ 2.
Indeed, a similar and with respect to the susceptibility even more significant effect is found
for the continuum states. In Figures 5.11 (population) and 5.12 (susceptibility) we present the distributions for states with small (panels (a): l ≤ 1) and large (panels (b): l ≥ 2) angular quantum
numbers. As noticed before for the excited states, it is the population in the states with larger
angular momentum quantum numbers which exhibit the nonadiabatic and nonperturbative characteristic features both in the population and in the contribution to the susceptibility. Interestingly,
in the case of the continuum the contributions from the higher l states are by far the dominant ones
at high intensities, giving rise to the signatures observed in the previous results (see Figure 5.7 forthe continuum-continuum contributions and Figure 5.5 for the total susceptibility).
The continuum population at the zeros of the field exhibits classic non-perturbative behavior
at 3 × 1014 W/cm2 : the growth in population is well described by an exponential curve. This
non-perturbative ionization corresponds with both the onset of a negative contribution to the

92
6

a)

b)

4
3

4
3

2

2

1

1

population (x10−5)

population (x10−6)

5

0

0
−8

−6

−4

−2

0
cycles

2

4

6

8

−8

−6

−4

−2

0
cycles24

68

Figure 5.11: Continuum population, for (a) l ≤ 1 and (b) l > 1 on the right. The lines represent results of calculations for 100 TWcm−2 (purple), 150 TWcm−2 (red), 200 TWcm−2 (green),
250 TWcm−2 (orange), and 300 TWcm−2 (blue) pulses. Note the difference in the time dependence
on the left, and on the rightsusceptibility and the breakdown in the susceptibility as shown in the previous Chapter, where at
3 × 1014 W/cm2 the ab-initio susceptibility separates from the perturbative susceptibility.

5.4

Conclusion
In this Chapter, we have extended our analysis of the nonlinear susceptibility towards its timedependence during the interaction of helium atom with a laser pulse. In our study, we focus on theintensity regime between 1 × 1014 and 3 × 1014 W/cm2 , in which we have previously identified thetransition from perturbative to nonperturbative light-matter interaction. In agreement with recent
experimental data, our theoretical results show several deviations from the adiabatic dependenceof the susceptibility on the field envelope that would be expected for a perturbative interaction.
In particular, we observe a shift of the maximum in the susceptibility towards times before thecenter of the pulse and a change in the slope of the results in the trailing part. We identify thatboth characteristic features are due to the onset of nonperturbative multiphoton ionization, thestrong population of the continuum and the transition of the population back from the continuum
into excited states. Moreover, our analysis shows that states with a higher angular quantum

93
number (l ≥ 1) are significantly populated at higher intensities, which provides a signature forthe nonperturbative interaction leading to the deviations from the typical adiabatic shape of thesusceptibility as a function of time.

94

a)

0.0

0.810

−0.2

0.805

−0.4

0.800

−0.6
−0.8

0.795

−1.0

0.790
0.785

−1.2

b)
−8

−6

−4

−2

0
cycles

2

4

6

8

−8

susceptibility (x10−2 a.u.)susceptibility (a.u.)

0.815

−6

−4

−2

0
cycles24

68

Figure 5.12: The total susceptibility, split up into l ≤ 1 on the left and l > 1 on the right. From
top to bottom, the lines represent 100TWcm−2 (purple), 150TWcm−2 (red), 200TWcm−2 (green),
250TWcm−2 (orange), and 300TWcm−2 (blue) pulses. The features we previously associated withthe continuum are completely absent in the l ≤ 1 susceptibility, instead they show up in the higher
l states.

Chapter 6

Conclusions

In this thesis, we have examined the response of matter to a strong light pulse, via an analysisof the electrical susceptibility of atoms at intensities in the transition regime between the perturbative and non-perturbative atom-field interaction. In Chapter 1 we discussed the introduction of thehigher order Kerr effect as a mechanism to explain the physics of filamentation. Filamentation is aphenomenon that occurs when propagating intense, short laser pulses in media. When the powerof the pulse is large enough, the beam will undergo self-focusing where the higher intensity parts ofthe beam will induce a higher index of refraction than the lower intensity parts. This self-focusing
will continue until ionization occurs, at which point the ionized electrons will defocus the beam
until the intensity has decreased enough to let self-focusing dominate again. The ionization anddefocusing of this process involve what have traditionally been understood as non-perturbative
processes – e.g. above threshold ionization, high harmonic generation and multiphoton ionization.
The process of filamentation has been long understood in terms of this focus, defocus cycle until
Loriot et al. [24] extracted out the terms of the index of refraction power series and found both
negative and positive terms. This led to a series of articles that discussed the implications of sucha series, including the possibility of ionization-free defocusing as a model of filamentation. This
discussion provides the motivation of this thesis, as it raised a more fundamental question regarding
the applicability of the perturbative series to a process that occurs simultaneously with well known
non-perturbative physics.
Examining the perturbative series description and non-perturbative version of the suscepti-

96bility side-by-side requires a detailed understanding of the perturbative series. In Chapter 2 wegave brief review of perturbation theory, a discussion of the first few perturbative dipole moment
terms and their corresponding susceptibilities, and some remarks about the relationship between
the susceptibility and the index of refraction. we discussed the relationship between the terms of
the perturbative index of refraction and the perturbative susceptibility, describing the nonlinear
dependence of the index of refraction on the density of the media, while the susceptibility’s linear
dependence on the density makes for a more appropriate quantity for single atom calculations. The
first few dipole moment terms describe linear polarization, sum and difference frequency generation,
along with the non-resonant physics that we are interested in. Finally, we discussed how to find thesusceptibilities from the dipole moment before describing the general form of the susceptibilities
for an arbitrary order in the perturbation series.
The perturbative susceptibilities require the field-free basis states for calculation, the methods of obtaining this basis comprised the bulk of Chapter 3. This Chapter is where we outlined our
method for finding the field-free energy basis states that we use for our ab-initio calculations alongwith our perturbative susceptibility calculations. The field-free basis was observed within a number of approximations to ease the numerical difficulties of solving the time dependent Schödinger
equation in its most general form. We applied a single-active electron approximation as well as
restricting ourselves to spherically symmetric potentials, in order to restrict ourselves to spherically separable problems in which we only need one dimension to solve numerically. From here,
we used the Numerov method to find the radial eigenstates, relying on the spherical harmonics forour angular states. The dimensionality of our problem was further reduced by limiting ourselves
to linear polarization of the external field. For the ab-initio method, we used a Crank-Nicholson
method for time propagation. This allowed us to find both the perturbative susceptibilities – using
the field-free basis and the general form of the susceptibilities in Chapter 2, and non-perturbative
ab-initio susceptibilities from the dipole moment using the same theoretical framework.
In Chapter 4, we made a comparison of the perturbative susceptibilities and the non-perturbative
ab-initio susceptibilities found using the ab-initio method. The perturbative series was shown to

97
be positive for all terms calculated for helium and hydrogen (up to the 11th order perturbative
correction or the sixth non-zero term in the series), which immediately countered the behavior seenin the higher order Kerr terms that Loriot et al. [24] found. We performed a comparison between
the perturbative series and the susceptibilities found using the non-perturbative ab-initio method
for a number of different intensities, and showed a clear threshold intensity for helium and hydrogen
where the perturbative and the non-perturbative susceptibilities diverge. At this intensity, we canclearly state that the perturbative series no longer describes the physics of the system. We next
extended this comparison to the third and fifth harmonics generated in hydrogen at longer wavelengths, where we found that the ab-initio calculation appeared to diverge from the perturbative
prediction at lower intensities for higher harmonics.
Finally, we discussed the time dependent susceptibility in Chapter 5. Here, we used theshort time Fourier transform to model the time dependent susceptibility. We examined the susceptibility of helium around the intensity that led to the divergence between the perturbative andnon-perturbative susceptibilities in the previous Chapter. Our analysis focused on the features ofthe susceptibility that drive this deviation. We split the susceptibility into continuum-continuum,
bound-bound and bound-continuum contributions using the field-free eigenstates to make the distinction. We showed that multiphoton ionization and the resultant ionized electrons in the continuum have the largest non-perturbative effect and cause the continuum-continuum susceptibility to
be negative. The bound-bound contribution is shown to be dominated by perturbative behavior.
This distinction was made more explicit by a separation into low (l < 2) and high (l ≥ 2) angularmomentum states. This separation showed that the perturbative behavior is confined to the low
angular momentum states, even when the continuum is taken into account, while the susceptibility from the high angular momentum states is primarily non-perturbative. The population of thecontinuum, split along these same lines, supported this idea.

Bibliography[1] T. H. Maiman. Stimulated Optical Radiation in Ruby. Nature, 187(4736):493–494, 1960.
[2] P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich. Generation of optical harmonics.
Physical Review Letters, 7(4):118–119, 1961.
[3] P. Agostini, F. Fabre, G. Mainfray, G. Petite, and N. K. Rahman. Free-free transitions following
six-photon ionization of xenon atoms. Physical Review Letters, 42(17):1127–1130, 1979.
[4] M. Hercher. Laser-induced damage in transparent media. Journal of the Optical Society ofAmerica, 54(563), 1964.
[5] G. Steinberg. Filamentary Tracks Formed in Transparent Optical Glass by Laser Beam SelfFocusing. I. Experimental Investigation*. Physical Review A, 4(3):1182–1194, 1971.
[6] E. L. Kerr. Filamentary Tracks Formed in Transparent Optical Glass by Laser Beam SelfFocusing. II. Theoretical Analysis*. Physical Review A, 4(3):1195–1218, 1971.
[7] C. G. Morgan. Laser-induced breakdown of gases. Reports on Progress in Physics, 38(5):621–
665, 1975.
[8] A. Braun, G. Korn, X. Liu, D. Du, J. Squier, and G. Mourou. Self-channeling of high-peakpower femtosecond laser pulses in air. Optics Letters, 20(1):73–75, 1995.
[9] E. T. Nibbering, P. F. Curley, G. Grillon, B. S. Prade, M. A. Franco, F. Salin, and A. Mysyrowicz. Conical emission from self-guided femtosecond pulses in air. Optics Letters, 21(1):62–65,
1996.
[10] B. La Fontaine, F. Vidal, Z. Jiang, C. Y. Chien, D. Comtois, A. Desparois, T. W. Johnston, J.C. Kieffer, H. Pépin, H. P. Mercure, H. Pepin, and H. P. Mercure. Filamentation of ultrashortpulse laser beams resulting from their propagation over long distances in air. Physics ofPlasmas, 6(5):1615, 1999.
[11] G. Méchain, A. Couairon, Y.-B. André, C. D’Amico, M. Franco, B. Prade, S. Tzortzakis,
A. Mysyrowicz, and R. Sauerbrey. Long-range self-channeling of infrared laser pulses in air: a
new propagation regime without ionization. Applied Physics B, 79(3):379–382, 2004.
[12] J. E. Rothenberg. Space-time focusing: breakdown of the slowly varying envelope approximation in the self-focusing of femtosecond pulses. Optics Letters, 17(19):1340, 1992.

99[13] C. P. Hauri, W. Kornelis, F. W. Helbing, A. Heinrich, A. Couairon, A. Mysyrowicz, J. Biegert,and U. Keller. Generation of intense, carrier-envelope phase-locked few-cycle laser pulses
through filamentation. Applied Physics B: Lasers and Optics, 79(6):673–677, 2004.
[14] S. Tzortzakis, G. Méchain, G. Patalano, Y.-B. André, B. Prade, M. Franco, A. Mysyrowicz,
J.-M. Munier, M. Gheudin, G. Beaudin, and P. Encrenaz. Coherent subterahertz radiationfrom femtosecond infrared filaments in air. Optics Letters, 27(21):1944–1946, 2002.
[15] A. Couairon and A. Mysyrowicz. Femtosecond filamentation in transparent media. Physics
Reports, 441(2-4):47–189, 2007.
[16] A. Becker, N. Aközbek, K. Vijayalakshmi, E. Oral, C.M. Bowden, and S.L. Chin. Intensity
clamping and re-focusing of intense femtosecond laser pulses in nitrogen molecular gas. Applied
Physics B, 73(3):287–290, 2001.
[17] M. Ferray, A. L’Huillier, X. F. Li, L. A. Lompre, G. Mainfray, and C. Manus. Multipleharmonic conversion of 1064 nm radiation in rare gases. Journal of Physics B: Atomic,
Molecular and Optical Physics, 21(3):L31–L35, 1999.
[18] R. W. Boyd. Nonlinear Optics, Third Edition. Academic Press, 3rd edition, 2008.
[19] J. H. Marburger. Self-focusing: Theory. Progress in Quantum Electronics, 4:35–110, 1975.
[20] Á Börzsönyi, Z. Heiner, A. P. Kovács, M. P. Kalashnikov, and K. Osvay. Measurement ofpressure dependent nonlinear refractive index of inert gases. Optics Express, 18(25):25847–
25854, 2010.
[21] H. Lange, A. Chiron, J.-F. Ripoche, A. Mysyrowicz, P. Breger, and P. Agostini. High-Order
Harmonic Generation and Quasiphase Matching in Xenon Using Self-Guided Femtosecond
Pulses. Physical Review Letters, 81(8):1611–1613, 1998.
[22] P. B. Corkum. Plasma perspective on strong field multiphoton ionization. Physical Review
Letters, 71(13):1994–1997, 1993.
[23] F. H. M. Faisal. Theory of Multiphoton Processes. Springer Science+Business Media, LLC,
New York, New York, USA, 1987.
[24] V. Loriot, E. Hertz, O. Faucher, and B. Lavorel. Measurement of high order Kerr refractive
index of major air components. Optics Express, 17(16):13429–13434, 2009.
[25] P. Béjot, J. Kasparian, S. Henin, V. Loriot, T. Vieillard, E. Hertz, O. Faucher, B. Lavorel, andJ. P. Wolf. Higher-Order Kerr Terms Allow Ionization-Free Filamentation in Gases. PhysicalReview Letters, 104(10):103903, 2010.
[26] P. Béjot, E. Hertz, J. Kasparian, B. Lavorel, J. P. Wolf, and O. Faucher. Transition from
Plasma-Driven to Kerr-Driven Laser Filamentation. Physical Review Letters, 106(24):243902,
2011.
[27] W. Ettoumi, P. Béjot, Y. Petit, V. Loriot, E. Hertz, O. Faucher, B. Lavorel, J. Kasparian, andJ. P. Wolf. Spectral dependence of purely-Kerr-driven filamentation in air and argon. PhysicalReview A, 82(3):33826, 2010.

100[28] V. Loriot, P. Béjot, W. Ettoumi, Y. Petit, J. Kasparian, S. Henin, E. Hertz, B. Lavorel,
O. Faucher, and J. P. Wolf. On negative higher-order Kerr effect and filamentation. Laser
Physics, 21(7):1319–1328, 2011.
[29] C. Köhler, R. Guichard, E. Lorin, S. Chelkowski, A. D. Bandrauk, L. Bergé, and S. Skupin.
Saturation of the nonlinear refractive index in atomic gases. Physical Review A, 87(4):43811,
2013.
[30] E. A. Volkova, A. M. Popov, and O. V. Tikhonova. Nonlinear polarization response of an
atomic gas medium in the field of a high-intensity femtosecond laser pulse. JETP Letters,
94(7):519–524, 2011.
[31] C. Bree, A. Demircan, and G. Steinmeyer. Method for Computing the Nonlinear Refractive
Index via Keldysh Theory. IEEE Journal of Quantum Electronics, 46(4):433–437, 2010.
[32] C. Bree, A. Demircan, and G. Steinmeyer. Saturation of the All-Optical Kerr Effect. PhysicalReview Letters, 106(18):183902, 2011.
[33] M. Kolesik, E. M. Wright, and J. V. Moloney. Femtosecond filamentation in air and higherorder nonlinearities. Optics Letters, 35(15):2550–2552, 2010.
[34] J. Ni, J. Yao, B. Zeng, W. Chu, G. Li, H. Zhang, C. Jing, S. L. Chin, Y. Cheng, and Z. Xu.
Comparative investigation of third- and fifth-harmonic generation in atomic and moleculargases driven by midinfrared ultrafast laser pulses. Physical Review A, 84(6):063846, 2011.
[35] G. O. Ariunbold, P. Polynkin, and J. V. Moloney. Third and fifth harmonic generation by
tightly focused femtosecond pulses at 2.2 µm wavelength in air. Optics Express, 20(2):1662–
1667, 2012.
[36] P. N. Butcher and D. Cotter. The Elements of Nonlinear Optics. Cambridge University Press,
New York, New York, USA, 1990.
[37] K. C. Kulander, F. H. Mies, and K. J. Schafer. Model for studies of laser-induced nonlinear
processes in molecules. Physical Review A, 53(4):2562–2570, 1996.
[38] M. P. de Boer and H. G. Muller. Observation of large populations in excited states after
short-pulse multiphoton ionization. Physical Review Letters, 68(18):2747–2750, 1992.
[39] I. A. Ivanov and A. S. Kheifets. Harmonic generation for atoms in fields of varying ellipticity:
Single-active-electron model with Hartree-Fock potential. Physical Review A, 79(5):053827,
2009.
[40] J. B. Krieger, Y. Li, and G. J. Iafrate. Systematic approximations to the optimized effective
potential: Application to orbital-density-functional theory. Physical Review A, 46(9):5453–
5458, 1992.
[41] X.-M. Tong and S.-I. Chu. Density-functional theory with optimized effective potential andself-interaction correction for ground states and autoionizing resonances. Physical Review A,
55(5):3406–3416, 1997.
[42] M. Abu-samha and L. B. Madsen. Single-active-electron potentials for molecules in intenselaser fields. Physical Review A, 81(3):033416, 2010.

101[43] X. M. Tong and C. D. Lin. Empirical formula for static field ionization rates of atoms andmolecules by lasers in the barrier-suppression regime. Journal of Physics B: Atomic, Molecularand Optical Physics, 38(15):2593–2600, 2005.
[44] P. M. Morse and H. Feshbach. Methods of theoretical physics. International series in pure andapplied physics. McGraw-Hill, New York, NY, 1953.
[45] T. E. Simos. A Numerov-type method for the numerical solution of the radial Schrödinger
equation. Applied Numerical Mathematics, 7(2):201–206, 1991.
[46] G. Beylkin, M. J. Mohlenkamp, and F. Pérez. Approximating a wavefunction as an unconstrained sum of Slater determinants. Journal of Mathematical Physics, 49(3):032107, 2008.
[47] C. J. Cramer. Essentials of Computational Chemistry: Theories and Models. Wiley, 2013.
[48] R. J. LeVeque. Finite Difference Methods for Ordinary and Partial Differential Equations.
Steady-State and Time-Dependent Problems. SIAM, 2007.
[49] L. Jianchun, G. A. Pope, and K. Sepehrnoori. A high-resolution finite-difference scheme fornonuniform grids. Applied mathematical modelling, 19(3):162–172, 1995.
[50] E. Tadmor. A Review of Numerical Methods for Nonlinear Partial Differential Equations.
Bulletin of the American Mathematical Society, 49(4):507–554, 2012.
[51] M. Soriano and J. J. Palacios. Theory of projections with non-orthogonal basis sets: Partitioning techniques and effective Hamiltonians. Physical Review B, 90(7):075128, 2014.
[52] H. A. Bethe and E. E. Salpeter. Quantum mechanics of one- and two-electron atoms. Plenum,
New York, NY, 1977.
[53] H. Bachau, E. Cormier, P. Decleva, J. E. Hansen, and F. Martı́n. Applications of B-splines inatomic and molecular physics. Reports on Progress in Physics, 64(12):1815–1943, 2001.
[54] C. de Boor. A Practical Guide to Splines. Springer, 2001.
[55] F. Martı́n. Ionization and dissociation using B-splines: Photoionization of the hydrogenmolecule. Journal of Physics B: Atomic, Molecular and Optical Physics, 32(16):R197–R231,
1999.
[56] O. Zatsarinny and K. Bartschat. The B-spline R-matrix method for atomic processes: application to atomic structure, electron collisions and photoionization. Journal of Physics B:
Atomic, Molecular and Optical Physics, 46(11):112001, 2013.
[57] J. C. Slater. Atomic Shielding Constants. Physical Review, 36(1):57, 1930.
[58] S. Datta. Evaluation of Coulomb integrals with hydrogenic and Slater-type orbitals. Journalof Physics B: Atomic and Molecular Physics, 18(5):853–857, 1985.
[59] W. J. Hehre, R. F. Stewart, and J. A. Pople. Self-Consistent Molecular-Orbital Methods. I.
Use of Gaussian Expansions of Slater-Type Atomic Orbitals. The Journal of Chemical Physics,
2003.

102
[60] V. P. Majety, A. Zielinski, and A. Scrinzi. Photoionization of few electron systems: a hybrid
coupled channels approach. New Journal of Physics, 17(6):063002, 2015.
[61] R. Kopold, W. Becker, M. Kleber, and G. G. Paulus. Channel-closing effects in high-order
above-threshold ionization and high-order harmonic generation. Journal of Physics B: Atomic,
Molecular and Optical Physics, 35(2):217–232, 2002.
[62] B. Yang, K. J. Schafer, B. Walker, K. C. Kulander, L. F. DiMauro, and P. Agostini. UP, 3UP,
11UP: Above-Threshold Ionization Revisited. Acta Physica Polonica A, 86(1-2):41–50, 1994.
[63] V. Pupyshev. The Lenz vector in the confined hydrogen atom problem. Chemical Physics
Letters, 295(3):217–222, 1998.
[64] I. A. Stegun and M. Abramowitz. Handbook of mathematical functions with formulas, graphs,
and mathematical tables. Dover books on intermediate advanced mathematics. U.S. National
Bureau of Standards, New York, NY, 1964.
[65] M. Galassi, B. Gough, F. Rossi, M. Booth, G. Jungman, J. Theiler, and J. Davies.
GNU Scientific Library. Reference Manual. Network Theory., 2001.
[66] L. N. Trefethen and D. Bau. Numerical Linear Algebra. Society for Industrial and Applied
Mathematics, Philadelphia, PA, 1997.
[67] N. Michel, Bogoliubov Laboratory of Theoretical Physics JINR, J.-C. Str, Dubna, and Russia.
Precise Coulomb wave functions for a wide range of complex `, η and z. Computer Physics
Communications, 176(3):232–249, 2007.
[68] A. Meurer, M. Paprocki, O. Čertı́k, M. Rocklin, S. Singh, R. P. Muller, F. Bonazzi, H. Gupta,
S. Vats, F. Johansson, F. Pedregosa, M. J. Curry, A. R. Terrel, Š. Roučka, A. Saboo, I. Fernando, S. Kulal, R. Cimrman, and A. Scopatz. Sympy: Symbolic computing in python. PeerJ
Computer Science, 3:e103, 2017.
[69] S. Chen, X. Gao, J. Li, A. Becker, and A. Jaroń-Becker. Application of a numerical-basisstate method to strong-field excitation and ionization of hydrogen atoms. Physical Review A,
86(1):013410–9, 2012.
[70] O. A. Sharafeddin, D. J. Kouri, and D. K. Hoffman. Time-dependent treatment of scattering:potential referenced and kinetic energy referenced modified Cayley approaches to atom–diatom
collisions and time-scale separations. Canadian Journal of Chemistry, 70(2):686–692, 2011.
[71] J. Crank and P. Nicolson. A practical method for numerical evaluation of solutions of partial
differential equations of the heat-conduction type [reprint of MR0019410 (8,409b)]. Advancesin Computational Mathematics, 6(3-4):207–226 (1997), 1996.
[72] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F. Smith,
S. Zampini, and H. Zhang. PETSc Web page. http://www.mcs.anl.gov/petsc, 2016.
[73] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F. Smith,
S. Zampini, and H. Zhang. PETSc users manual. Technical Report ANL-95/11 - Revision 3.7,
Argonne National Laboratory, 2016.

103[74] S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith. Efficient management of parallelism inobject oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen,editors, Modern Software Tools in Scientific Computing, pages 163–202. Birkhäuser Press,
1997.
[75] M. Nurhuda, A. Suda, and K. Midorikawa. Saturation of dynamic nonlinear susceptibility of
noble gas atom in intense laser field. RIKEN REVIEW, (48):40–43, 2002.
[76] M. Nurhuda, A. Suda, and K. Midorikawa. Generalization of the Kerr effect for high intensity,
ultrashort laser pulses. New Journal of Physics, 10(5):053006, 2008.
[77] A. McPherson, G. Gibson, H. Jara, U. Johann, T. S. Luk, I. A. McIntyre, K. Boyer, and C. K.
Rhodes. Studies of multiphoton production of vacuum-ultraviolet radiation in the rare gases.
Journal of the Optical Society of America B, 4(4):595–601, 1987.
[78] X. F. Li, A. L’Huillier, M. Ferray, L. A. Lompre, and G. Mainfray. Multiple-harmonic generation in rare gases at high laser intensity. Physical Review A, 39(11):5751–5761, 1989.
[79] A. Spott, A. Jaroń-Becker, and A. Becker. Ab initio and perturbative calculations of theelectric susceptibility of atomic hydrogen. Physical Review A, 90(1):013426, 2014.
[80] F. Eilenberger, M. Bache, and S. Minardi. Higher-order Kerr effect and harmonic cascadingin gases. Optics Letters, 37(22):4612–4614, 2012.
[81] J. M. Brown, E. M. Wright, J. V. Moloney, and M. Kolesik. On the relative roles of higher-ordernonlinearity and ionization in ultrafast light-matter interactions. Optics Letters, 37(10):1604–
1606, 2012.
[82] 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(11):815–820, 2009.
[83] C. H. Garcı́a. Coherent attosecond light sources based on high-order harmonic generation:
influence of the propagation effects. PhD thesis, Salamanca, 2013.
[84] A. Rundquist, C. G. Durfee, Z Chang, C. Herne, S. Backus, M. M. Murnane, and H. C.
Kapteyn. Phase-Matched Generation of Coherent Soft X-rays. Science, 280(5368):1412–1415,
1998.
[85] J. A. Armstrong, N. Bloembergen, J. Ducuing, and P. S. Pershan. Interactions between Light
Waves in a Nonlinear Dielectric. Physical Review, 127(6):1918–1939, 1962.
[86] C. Hernandez-Garcia, T. Popmintchev, M. M. Murnane, H. C. Kapteyn, L. Plaja, A. Becker,
and A. Jaroń-Becker. Group velocity matching in high-order harmonic generation driven by
mid-infrared lasers. New Journal of Physics, 18(7):073031–10, 2016.
[87] T. Popmintchev, M.-C. Chen, A. Bahabad, M. Gerrity, P. Sidorenko, O. Cohen, M. M. Christov, I. P. Murnane, and H. C. Kapteyn. Phase matching of high harmonic generation in thesoft and hard X-ray regions of the spectrum. Proceedings of the National Academy of Sciences,106(26):10516–10521, 2009.

104[88] 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(12):822–832, 2010.
[89] T. Popmintchev, M. C. Chen, D. Popmintchev, P. Arpin, S. Brown, S. Alisauskas, G. Andriukaitis, T. Balciunas, O. D. Mucke, A. Pugzlys, A. Baltuska, B. Shim, S. E. Schrauth,A. Gaeta, C. Hernandez-Garcia, L. Plaja, A. Becker, A. Jaroń-Becker, M. M. Murnane, and
H. C. Kapteyn. Bright Coherent Ultrahigh Harmonics in the keV X-ray Regime from MidInfrared Femtosecond Lasers. Science, 336(6086):1287–1291, 2012.
[90] D. Popmintchev, C. Hernandez-Garcia, F. Dollar, C. Mancuso, J. A. Perez-Hernandez, M. C.
Chen, A. Hankla, X. Gao, B. Shim, A. L. Gaeta, M. Tarazkar, D. A. Romanov, R. J. Levis,
J. A. Gaffney, M. Foord, S. B. Libby, A. Jaroń-Becker, A. Becker, L. Plaja, M. M. Murnane,
H. C. Kapteyn, and T. Popmintchev. Ultraviolet surprise: Efficient soft x-ray high-harmonic
generation in multiply ionized plasmas. Science, 350(6265):1225–1231, 2015.
[91] J. K. Wahlstrand, Y. H. Cheng, and H. M. Milchberg. High Field Optical Nonlinearity and
the Kramers-Kronig Relations. Physical Review Letters, 109(11):113904, 2012.
[92] J. K. Wahlstrand, Y. H. Cheng, Y. H. Chen, and H. M. Milchberg. Optical Nonlinearity in Ar
and N {2} near the Ionization Threshold. Physical Review Letters, 107(10):103901, 2011.
[93] J. K. Wahlstrand and H. M. Milchberg. Effect of a plasma grating on pump-probe experiments
near the ionization threshold in gases. Optics Letters, (3):3, 2011.
[94] J. Allen. Short term spectral analysis, synthesis, and modification by discrete Fourier transform. IEEE Transactions on Acoustics, Speech, and Signal Processing, 25(3):235–238, 1977.
[95] G. DAntona and A. Ferrero. Digital Signal Processing for Measurement Systems. Theory andApplications. Springer Science & Business Media, 2006.
[96] F. J. Harris. On the use of windows for harmonic analysis with the discrete Fourier transform.
Proceedings of the IEEE, 66(1):51–83, 1978.
[97] P. Lynch. The Dolph–Chebyshev window: A simple optimal filter. Monthly weather review,
125(4):655–660, 1997.
[98] R. R. Jones, D. W. Schumacher, and P. H. Bucksbaum. Population trapping in Kr and Xe in
intense laser fields. Physical Review A, 47(1):R49–R52, 1993.
[99] T. Nubbemeyer, K. Gorling, A. Saenz, U. Eichmann, and W. Sandner. Strong-Field Tunnelingwithout Ionization. Physical Review Letters, 101(23):233001, 2008.