Atoms, Molecules and Clusters in Intense Laser
FieldsbyZachary B. Walters
B.S., Harvey Mudd College, 2002

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
2009This thesis entitled:
Atoms, Molecules and Clusters in Intense Laser Fields
written by Zachary B. Walters
has been approved for the Department of Physics

Chris H. Greene

Murray Holland

Date

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

iii
Walters, Zachary B. (Ph.D., Physics)
Atoms, Molecules and Clusters in Intense Laser Fields
Thesis directed by Chris H. Greene

Recent advances in the technology of intense, short laser pulses have opened thepossibility of investigating processes in atoms, molecules and clusters in which the normal intramolecular forces between electrons and nuclei, and between different electrons,
are rivaled in strength by interactions with the driving laser, or with a cluster plasma.
Experiments using rescattered electrons offer a means of probing atomic and molecular
processes on ultrafast timescales.
This thesis extends techniques and concepts of atomic and molecular physics todescribe physics in the strong field regime. This involves investigating how electron
scattering from atoms and molecules is affected by the intense and time-varying electricfield of the laser, the effect of such scattering on experimental observables, and the roleof intramolecular structure on strong field processes.
Also investigated is the evolution of van der Waals atomic clusters when subject tointense laser pulses in the VUV regime. Here processes such as photoionization, inverse
bremsstrahlung heating, and collisional ionization and recombination are affected bothby the non-hydrogenic nature of the relevant atomic potentials but also by the screeningof these potentials by the cluster plasma.

Dedication

To my family.

v

Acknowledgements

Many people helped make this thesis possible, and I would like to express mygratitude to them all. More than anyone, I must thank my advisor, Chris Greene, formany patient years as a teacher and mentor. Along with his unparalleled knowledge ofatomic and molecular physics, Chris has a wonderful grasp of the subtle art of turninga physical insight into a physical theory.
Physics can be a talkative subject, and there are many people I would like to thankon this front. Robin Santra, Seth Rittenhouse, Nirav Mehta, Jose D’Incao, Stefano
Tonzani, Roman Curik, Dan Haxton, Javier von Stecher, Jia Wang, Josh Dunn and
Heather Crowell were all members of the Greene group during my graduate career, andI would like to thank each of them for many animated hours in front of a whiteboard
during my graduate career. Ed Meyer also deserves thanks for his helpful explanationsof quantum chemistry.
Much of the work in this thesis was inspired by experiments performed by thegroup of Henry Kapteyn and Margaret Murnane. Aside from Henry and Margaret, I
would like to thank Nick Wagner, Andrea Wuest, Robynne Lock, Wen Li and Xibin
Zhou for many stimulating and informative discussions.
Pete Ruprecht had the sometimes thankless task of keeping all the computers
working, and has been a fount of useful information about programming issues.
Finally I would like to thank Pam Leland for her administrative support.

vi
The funding for this research was provided by the Department of Energy and the
National Science Foundation Extreme Ultraviolet Engineering Research Center. Computational support was provided by the W.M. Keck Foundation.

vii

Contents

Chapter
1 Introduction and Physical Principles

1

2 Scattering Resonances in the Presence of Intense Laser Fields

5

2.1

Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . .

18

3 Molecular scattering and the Limits of the Plane Wave Approximation in the
Measurement of Molecular Properties

21

3.1

Scattering States and Ramifications for Molecular Tomography . . . . .

22

3.2

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

3.3

Appendix: Gauges and Dispersion Relations . . . . . . . . . . . . . . . .

40

4 High Harmonic Generation in SF6 : Raman-excited Vibrational Quantum Beats 43
4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

4.2

The Vibrational Wavefunction of the Molecule . . . . . . . . . . . . . .

44

4.3

Vibrational Interference . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

4.4

Describing the Continuum Electron . . . . . . . . . . . . . . . . . . . . .

52

4.4.1

Tunneling Ionization . . . . . . . . . . . . . . . . . . . . . . . . .

54

4.4.2

Semiclassical Propagation . . . . . . . . . . . . . . . . . . . . . .

56

4.4.3

Stationary phase calculation of scattering coefficients . . . . . . .

58

Comparison with Experiment . . . . . . . . . . . . . . . . . . . . . . . .

62

4.5

viii
4.6

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

4.7

Appendix: Tunneling Ionization . . . . . . . . . . . . . . . . . . . . . . .

65

5 Laser-cluster interactions in the VUV energy range

67

5.1

Photoionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

5.2

Inverse Bremsstrahlung Heating . . . . . . . . . . . . . . . . . . . . . . .

72

5.3

Collisional Ionization and Recombination . . . . . . . . . . . . . . . . .

76

5.4

Cluster Dynamics during the Laser Pulse . . . . . . . . . . . . . . . . .

83

5.5

Nonideal Plasma Screening . . . . . . . . . . . . . . . . . . . . . . . . .

89

5.6

Hydrogenic Model of Inverse Bremsstrahlung . . . . . . . . . . . . . . .

97

5.7

Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . 102

6 Summary, Conclusions and Outlook

106

Bibliography

108

ix

FiguresFigure2.1

Adiabatic potential energy curves for Cs+ , calculated for E0 = 0 and
E0 = .02. The potential curves depart from each other at long range, but
are largely unaffected in the region of the f-wave potential energy well. .

2.2

Comparison between recombination dipole calculated in this work andby Krause et al [28]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.3

7

10

In the vicinity of xenon, the f wave potential falls rapidly with increasing
nuclear charge. At lanthanum (Z=57), the potential supports a boundstate. Here f wave Herman-Skillman potentials, are plotted for a variety
of xenon-like atoms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.4

12

a) Photoionization cross section vs energy for the 5d subshell of Cs + in
the 5p5 5d1 excited state. b) p- and f-wave scattering phase shifts vsenergy. The shape resonance corresponds with a change in the f-wave
phase shift of π. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.5

13

a) Photoionization cross section vs energy for the 6s subshell of Cs + in
the 5p5 6s1 excited state. b) p-wave scattering phase shifts vs energy for
Cs+ . In contrast to the 5d shell, no shape resonance effects are present.

14

x
2.6

a) Cs+ 5d single channel photoionization cross section vs energy at low
electric field strengths. b)Cs+ 6s single channel photoionization crosssection vs energy at low electric field strengths. c) ratio of 5d to 6s
photoionization cross sections vs energy at low electric field strengths. .

2.7

16

a) Cs+ 5d single channel photoionization cross section vs energy at fieldstrengths comparable to pulse peak b)Cs+ 6s single channel photoionization cross section vs energy at field strengths comparable to pulse peakstrengths. c) ratio of 5d to 6s photoionization cross sections vs energy at
field strengths comparable to pulse peak.

2.8. . . . . . . . . . . . . . . . .

17a) ratio of interpolated d harmonic spectrum to interpolated s harmoniccross section, normalized by ADK ionization rate b) difference in phasebetween d and s harmonic spectra. . . . . . . . . . . . . . . . . . . . . .

3.1

19

The tomographic reconstruction procedure applied to the 1D square well.
(Top) Comparison of dipole matrix elements dk = hψc (x)| x |ψg (x)i, calculated using plane waves and scattering states for ψc (x). (Bottom) Because the scattering state matrix elements differ from those calculatedusing plane waves, the reconstructed image of the orbital will differ fromthe true bound state wavefunction. . . . . . . . . . . . . . . . . . . . . .

3.2

26

N2 photoionization cross sections vs photon energy. Calculations madeusing FERM3D and the plane wave approximation are compared to experimental measurements taken from [6]. . . . . . . . . . . . . . . . . . .

3.3

28

F2 photoionization cross sections vs photon energy. Calculations madeusing FERM3D and the plane wave approximation are compared to theoretical calculations taken from [27]. . . . . . . . . . . . . . . . . . . . .

28

xi
3.4

Comparison of the Hartree-Fock orbital and associated tomographic images for the N2 1πu orbital. a) and b) give the real and (zero) imaginary
components of the Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized
dipole matrix element. e) and f) give the real and imaginary components of the tomographic image made from the y-polarized dipole matrix
element. g) and h) give the real and imaginary components of the tomographic image made from the z-polarized dipole matrix element. . . . .

3.5

31

Comparison of the Hartree-Fock orbital and associated tomographic images for the N2 3σg orbital. a) and b) give the real and (zero) imaginary
components of the Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized
dipole matrix element. e) and f) give the real and imaginary components of the tomographic image made from the y-polarized dipole matrix
element. g) and h) give the real and imaginary components of the tomographic image made from the z-polarized dipole matrix element. . . . .

3.6

32

Comparison of the Hartree-Fock orbital and associated tomographic images for the N2 2σu orbital. a) and b) give the real and (zero) imaginary
components of the Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized
dipole matrix element. e) and f) give the real and imaginary components of the tomographic image made from the y-polarized dipole matrix
element. g) and h) give the real and imaginary components of the tomographic image made from the z-polarized dipole matrix element. . . . .

33

xii
3.7

Comparison of the Hartree-Fock orbital and associated tomographic images for the N2 2σg orbital. a) and b) give the real and (zero) imaginary
components of the Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized
dipole matrix element. e) and f) give the real and imaginary components of the tomographic image made from the y-polarized dipole matrix
element. g) and h) give the real and imaginary components of the tomographic image made from the z-polarized dipole matrix element. . . . .

3.8

34

Comparison of the Hartree-Fock orbital and associated tomographic images for the F2 1πg orbital. a) and b) give the real and (zero) imaginary
components of the Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized
dipole matrix element. e) and f) give the real and imaginary components of the tomographic image made from the y-polarized dipole matrix
element. g) and h) give the real and imaginary components of the tomographic image made from the z-polarized dipole matrix element. . . . .

3.9

35

Comparison of the Hartree-Fock orbital and associated tomographic images for the F2 3σg orbital. a) and b) give the real and (zero) imaginary
components of the Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized
dipole matrix element. e) and f) give the real and imaginary components of the tomographic image made from the y-polarized dipole matrix
element. g) and h) give the real and imaginary components of the tomographic image made from the z-polarized dipole matrix element. . . . .

36

xiii
3.10 Comparison of the Hartree-Fock orbital and associated tomographic images for the F2 1πu orbital. a) and b) give the real and (zero) imaginary
components of the Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized
dipole matrix element. e) and f) give the real and imaginary components of the tomographic image made from the y-polarized dipole matrix
element. g) and h) give the real and imaginary components of the tomographic image made from the z-polarized dipole matrix element. . . . .

37

3.11 Comparison of the Hartree-Fock orbital and associated tomographic images for the F2 2σu orbital. a) and b) give the real and (zero) imaginary
components of the Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized
dipole matrix element. e) and f) give the real and imaginary components of the tomographic image made from the y-polarized dipole matrix
element. g) and h) give the real and imaginary components of the tomographic image made from the z-polarized dipole matrix element. . . . .

38

3.12 Comparison of the Hartree-Fock orbital and associated tomographic images for the F2 2σg orbital. a) and b) give the real and (zero) imaginary
components of the Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized
dipole matrix element. e) and f) give the real and imaginary components of the tomographic image made from the y-polarized dipole matrix
element. g) and h) give the real and imaginary components of the tomographic image made from the z-polarized dipole matrix element. . . . .

39

xiv
3.13 a) Magnitude of the overlap between (normalized) tomographic images ofa bound wavefunction and the true wavefunction, and between differentptomographic images, calculated using q(²) = 2(k 2 /2 − ²V ), ² ∈ [0, 1].
b) Comparison of the maximally overlapping tomographic images to thetrue wavefunction. In the momentum gauge, maximal overlap was obtained for ² = .26, while in the length gauge, maximal overlap was obtained for ² = 1. Both images have been normalized and rotated to givea purely real overlap with the true wavefunction. . . . . . . . . . . . . .
4.1

42

The vibrational interference model [65] in one dimension. The moleculeends the first (Raman) pulse in a superposition of the ν = 0 and ν = 1
vibrational states. After a time delay, the two vibrational states are
mixed by stimulated Raman scattering (transfer matrix M), “hopping”
during ionization (I) and recombination (R), as well as evolution of theionic wavefunction while the electron is away (N). Interference between
adjacent vibrational states modulates the high harmonic signal. . . . . .

4.2

46

a)Spherically symmetric A1g (breathing mode) distortions change electronic state energies, but preserve the triple degeneracy. Nontotally symmetric Eg (b) and T2g (c) distortions break the triple degeneracy of SF+
6at the maximum symmetry point. Adiabatic energies are fit to the eigenvalues of the vibronic coupling matrix (equation 4.9) to solve for thevibronic coupling constants. . . . . . . . . . . . . . . . . . . . . . . . . .

4.3

51

One of three degenerate orbitals of SF6 . Red denotes positive lobes; blue
denotes negative lobes. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

xv
4.4

The tunneling wavefunction is approximated as an unperturbed molecular HOMO inside the classically allowed region, connecting to a declining
WKB exponential in the classically forbidden region. Stationary phasetrajectories leave from the outer turning point, beginning with zero velocity at time of ionization. . . . . . . . . . . . . . . . . . . . . . . . . .

4.5

56

Peak-to-peak modulation of the high harmonic signal vs. wavenumber,comparing theory to the two experimental runs for which data is available.
For the bars labeled “phase information included”, the angular average
was performed allowing δ in equation 4.17 to vary as a function of angle.
For the bars labeled “phase information excluded,” δ was set to zero for
all angles. Modulations corresponding to the same frequency have beenplaced side-to-side for purpose of comparison. (Figure taken from [65]) .

5.1

64

Inverse bremsstrahlung cross sections [Eqs. 5.13 and 5.14] for an electronwith incident energy E to absorb a 12.7-eV photon are given for an electron in the field of a purely Coulombic 1+ potential and for an electronin the field of a Xe Herman-Skillman atomic potential. The effects ofatomic structure on inverse bremsstrahlung rates are quite pronounced.
Figure taken from [64]. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.2

77

The inverse bremsstrahlung cross section as a function of energy is shown
for an electron in the field of a Debye-screened Xe Herman-Skillman potential, with the Debye screening length λD ranging from 1 a.u. to 20
a.u. As λD grows, the cross section approaches the limit of no plasma
screening, shown in this graph by the dotted line. As the Debye lengthof the cluster plasma shrinks, the charged ion is shielded more effectively
from the scattering electron, and the inverse bremsstrahlung cross section
is decreased. Figure taken from [64]. . . . . . . . . . . . . . . . . . . . .

77

xvi
5.3

Evolution of a 1500 atom cluster exposed to a 100 fs, 7 × 10 13 W/cm2
pulse, employing only photoionization and inverse bremsstrahlung heating. a) Energy absorbed vs. time. b) Ionic population vs. time. Xe2+
and Xe3+ are produced efficiently via photoionization. Figure taken from[64]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.4

78

Comparison of thermalization time with time to undergo stimulated orinverse bremsstrahlung for an electron in the cluster plasma, calculated
for a 1500 atom cluster exposed to a 7 × 1013 W/cm2 , 100 fs pulse, using
the Wigner-Seitz cutoff model for screening. The electrons thermalize
much faster than they interact with the laser at all times during the laser
pulse. Figure taken from [64]. . . . . . . . . . . . . . . . . . . . . . . . .

5.5Ejected electron spectrum.

85

Comparison between data from [30] andspatially-averaged spectra calculated using 70 atom clusters exposed to a
4.4 × 1012 W/cm2 , 100 fs pulse for two different models of plasma screening. The Wigner Seitz cutoff model uses the ordinary Debye length as
the screening radius, but the screening radius is not allowed to fall below xenon’s Wigner-Seitz radius at liquid density, 4.64 bohr. The Attard
model of screening calculates the screening radius according to equation(5.34), discussed in Section 5.5. The spectrum calculated using xenon’s
Wigner-Seitz radius as a minimum screening distance displays a strong
similarity to the experimental curve. The intensity of the experimental
spectra is arbitrary; magnitudes were chosen by setting each curve equal
at the beginning of the exponential tail in the experiment. Figure takenfrom [64]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90

xvii
5.6

Ejected electron spectra, calculated for the two sets of parameters andthe two models of screening. The Wigner Seitz cutoff model uses theordinary Debye length as the screening radius, but the screening radius
is not allowed to fall below xenon’s Wigner-Seitz radius at liquid density,
4.64 bohr. The Attard model of screening calculates the screening radiusaccording to equation (5.34), discussed in section 5.5. a) Nature parameters: 1500 atom clusters exposed to a 100 fs, 7 × 1013 W/cm2 pulse.
b) Thesis parameters: 2500 atom clusters exposed to a 50 fs, 2.5 × 10 13
W/cm2 pulse. Since electrons faster than about 1 eV are ejected fromthe cluster during the pulse the ejection spectra could serve as a window
into the dynamics of the laser-cluster interaction. Figure taken from [64].

5.7

91

The effects of collisional ionization and recombination are to allow theformation of charge states beyond Xe3+ Pictured is the time evolution ofa single 1500 atom cluster exposed to a 100 fs, 7 × 1013 W/cm2 pulse.
These states enhance the rate of inverse bremsstrahlung heating. As
the plasma expands and cools, the chemical equilibrium shifts toward
lower charge states on a timescale much longer than the laser pulse, until
decreasing plasma density causes recombination and ionization rates to
go to zero. a) Energy absorbed vs. time b) Ionic population vs. time
during laser pulse. Figure taken from [64]. . . . . . . . . . . . . . . . . .

5.8

92

The interaction of plasma screening with atomic potentials is unknown
as the screening length becomes very short. Here the screening length vstime is given for two simulations of a 1500 atom cluster exposed to a 100
fs, 7 × 1013 W/cm2 pulse, using two models for screening. In the firstmodel, the screening length is not allowed to fall below xenon’s WignerSeitz radius at liquid density. The second model for screening uses aformula given by Attard. Figure taken from [64]. . . . . . . . . . . . . .

94

xviii
5.9

Near the center of the pulse, calculating the evolution of the cluster using
Attard screening yields screening lengths shorter than the Wigner-Seitz
cutoff, allowing easier formation of high charge states than when theevolution is calculated for Wigner-Seitz screening, shown in Figure 5.7.
For a 1500 atom cluster exposed to a 7 × 1013 W/cm2 , 100 fs pulse:a) Energy absorbed vs time for the Attard screening model, b) Charge
species population vs time for the Attard screening model. Figure takenfrom [64]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95

5.10 Plasma coupling parameter vs. time for the two models of plasma screen3/2

1
ing. The coupling parameters are defined by Γee = aTand Γei = ZΓee3
)1/3
where the average distance between electrons a is given by a = ( 4πneand Z is the average charge of the ions. The plasma becomes very strongly
coupled early in the pulse, but the strength of the coupling decreases as
the plasma absorbs more energy in the course of the cluster heating. a)
Coupling parameters vs. time using Wigner-Seitz cutoff. b) Coupling
parameters vs. time for the Attard screening model. Figure taken from[64]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96

5.11 Inverse bremsstrahlung cross sections [Eqs. 5.13 and 5.14] calculated for
an electron in the field of a purely Coulombic 5+ potential and for anelectron in the field of a Xe Herman-Skillman atomic potential of thesame charge. In comparison with figure 5.1, it can be seen that at higher
charge states, the impact of atomic structure on inverse bremsstrahlungcross sections is decreased. Figure taken from [64]. . . . . . . . . . . . .

99

xix
5.12 Simulation of the laser-cluster interaction using a physical model takenfrom [52]. In this model, inverse bremsstrahlung cross sections are calculated using hydrogenic potentials and all high charge states are produced
via sequential photoionization. Collisional ionization and recombination
are not considered. a)Energy absorbed vs time. b) Charge state population vs time. c) Debye length vs time. Figure taken from [64]. . . . . . . 100
5.13 Evolution of the cluster calculated in the limit of long screening length,
λD = 20 bohr for a 1500 atom cluster exposed to a 100 fs, 7 × 10 13
W/cm2 pulse. Weak screening prevents ionization potential lowering due
to plasma effects, and precludes photoionization past Xe+ . The effect ofneglecting plasma screening effects is to reduce formation of high charge
states and to reduce the total energy absorbed by the cluster. The pulse
and cluster parameters are identical to those used in Figures 5.3, 5.7 and5.9. a) Energy absorbed vs. time b) Ionic population vs. time duringlaser pulse. Figure taken from [64]. . . . . . . . . . . . . . . . . . . . . . 101

Chapter 1

Introduction and Physical Principles

Under normal circumstances, an atom or a molecule passes through space like apassenger train speeding through the night. As a passenger on the train, interactions
with other passengers are direct and immediate. No two passengers may sit in the sameseat, and one may be attracted to some passengers and repelled from others. Some
passengers may be liable to move about the cars, while others may spend the entire
journey dozing in the same seat.
In contrast, the outside world is experienced at one remove. If one side of the train
has a particularly stunning view, you might turn to look at it while still maintaining
your distance from the crying baby in front of you and your proximity to the interesting
conversationalist in the back. Likewise, if a fellow passenger disembarks you mightsettle happily into his seat while caring nothing about which station he left at or why.
Most of the time, atomic and molecular physics works in the same way. The forces
between electrons and nuclei, or between different electrons in the same atom completely
dominate anything the outside world can bring to bear.
Much of the interest of intense laser-target interactions is that the interaction
with the outside world – in this case with the field imposed by the laser – can rival
the strength of the interactions within the original system. To return to the train
metaphor, suppose the glare from the rising sun made sitting on one side of the car
extremely irritating. Even though all passengers relate to one another in the same way,2
the seating arrangement might end up entirely different. If the train were following
a switchback, causing the glare to alternate which sides of the train it affected, theresulting arrangement of passengers might be extremely complex and time dependent.
In contrast to the molecule-centric perspective of AMO physics, strong field
physics might occasionally be accused of ignoring atoms and molecules completely.
Primarily concerned with the ways in which electrons interact with the strong electromagnetic fields encountered in an intense laser beam, frequently used theories often
treat atoms and molecules as an afterthought, a generic coulomb potential which altersa few results in non-essential ways. Two examples of this will be seen in later chapters. In semiclassical theories of high harmonic generation, encountered in Chapter 4
the Coulomb potential slightly alters the trajectories which contribute to a particularharmonic, but play little role in the intensity of that harmonic or the energy at which
the spectrum dies off. If the laser should happen to hit a cluster and create a nanoscale
plasma as in Chapter 5, the treatment of that plasma is not particularly different if theindividual atoms are xenon than if they were hydrogen.
A major goal of this thesis has been to unite these two perspectives: to do atomic
and molecular physics in a regime where external forces rival internal forces in strength,
or strong field physics in a way where the individual characteristics of particular atoms
and molecules may be accounted for.
Because this is such a broad topic, no thesis can hope to describe all cases ofinterest. In the pages that follow, different aspects of the laser-target interaction will
be dealt with in problem-specific ways. A common theme will be that the identity ofa particular atom or molecule makes itself important when continuum electrons scatterfrom that particle. Far from the particle, there is little to distinguish one ion of aparticular charge from another; it is only when an electron feels a short range potentialthat the identity of the particle becomes significant.
In the region close to the molecule, internal forces dominate. The electronic part3of a molecular wavefunction will consist of many core electrons which are contained
in deeply bound orbitals, plus a few electrons contained in less deeply bound valenceorbitals. Because these valence electrons are screened from the nuclear potentials by
the core electrons, they are much freer to rearrange themselves when subject to the fieldof a strong laser.
It is readily apparent that a full treatment of the dynamics of such a system willrapidly become intractable. Because every new particle in the system introduces three
new spatial degrees of freedom, a simple nitrogen molecule will already require 14*3=42
dimensions for the electronic wavefunctions alone! For this reason, most current theories
of strong field physics prefer to work with as few active electrons as possible, andusually only the electron in the highest occupied molecular orbital (HOMO) is allowed
to respond to the laser field.
This lone electron will then feel a strong but local potential due to interactions
with the parent ion, plus a weaker (typically on the order of a few tenths or hundredthsof an atomic unit) but long range electric field applied by the driving laser. Far from themolecule, the oscillating electric field becomes the dominant force felt by the electron.
Because the effects of the parent ion are felt only at short range, the interactionsof molecules with intense lasers share many characteristics regardless of the choice ofmolecule. These characteristics can be explained semiclassically using the three stepmodel introduced by Paul Corkum [7, 33]. In this model, electrons enter the continuum
by tunneling away from the parent ion, then follow classical trajectories in the oscillating laser field before recolliding with the parent ion. Because the electron is being
accelerated by the electric field of the laser , it can have a much greater energy at thetime of recollision than at the time of tunneling. If the electron then recollides withthe parent ion, it can emit a photon with an energy much higher than the frequency ofthe driving laser. Alternatively, the electron may scatter elastically from the ion and be
accelerated once again by the driving laser, or it can scatter inelastically and liberate

4
additional electrons from the parent ion.
The three step model is at its best in describing the evolution of the electron
while it is far away from the parent ion, and in describing gross features of the experimental observables. Because the tunneling and rescattering both occur at the locationof the parent ion, there are calculable maximum energies of the electron at the time ofrecollision, and after scattering from the parent ion. However, the semiclassical approximations at the heart of the three step model will be at their worst in describing theelectron in the vicinity of the parent ion, where the potential of the parent ion becomesparamount and the classical trajectories become very complicated. More sophisticatedtreatments of scattering are necessary to adequately describe this terminal phase.
Chapters 2 and 3 are concerned with details of the rescattering process. Chapter
2 addresses the question of how scattering features such as shape resonances are affectedby the intense, time dependent field of the laser. Chapter 3 addresses the question of
how electron-molecule scattering affects efforts such as molecular tomography, whichattempt to reconstruct molecular properties from experimental data.
Chapter 4 moves beyond the single active electron model to consider the evolution of the molecule’s vibrational wavefunction during the high harmonic process.
In addition, the three step Corkum model is extended to connect with more sophisticated short range models of electron-molecule scattering, thereby helping to unify themolecule-centric and strong field-centric approaches to HHG.
Chapter 5 deals with the problem of an atomic cluster in an intense beam of xrays. In contrast with earlier chapters, here only single photon processes are important.
The complications arise due to the dense and rapidly heating cluster plasma whichdevelops during the course of the pulse. This plasma acts to screen the atomic electrons,
lowering the binding energies and affecting cross sections for photoionization and inverse
bremsstrahlung heating.

Chapter 2

Scattering Resonances in the Presence of Intense Laser Fields

A central feature of many current efforts to probe molecular properties on ultrafast
timescales is the use of rescattered electrons as probes. Rescattering [7] occurs when
an electron tunnels free from a parent ion, accelerates in the field of the driving laser,
then recollides with the parent ion. Information about the molecule may then be foundby monitoring various scattering output channels. If the electron recollides with theparent ion, its high kinetic energy results in the emission of very high energy photons.
Alternatively, the electron may scatter elastically from the target or cause the target tobreak up into two or more fragments. Understanding the scattering process is thus anecessary step in interpreting rescattering experiments.
Aside from its experimental relevance, the rescattering process is interesting in its
own right, as an exotic variant of scattering. Because rescattering necessarily takes placein the presence of an intense, time varying laser field, familiar results from zero field,
time independent calculations may be altered in the context of strong field physics. This
chapter investigates the familiar concepts of shape resonances and potential barriers in
the context of a rescattering experiment.
Commonly encountered in scattering physics, a shape resonance consists of aquasi-bound state at positive energy. A potential barrier traps the scattering particle
in a local potential well, where it may survive for a long time before tunneling throughthe barrier to the continuum. The scattering wavefunction will often have a very large

6amplitude in the region behind the potential barrier over large energy ranges, withcorrespondingly large effects on experimental observables such as photoionization or-recombination cross sections [55].
The presence of a strong laser field complicates the static view of shape resonances. In the static limit, potential barriers may be wiped out by the strength of theapplied field at times during the laser cycle. Furthermore, in a time dependent treatment, nonadiabatic coupling between scattering channels could limit the buildup of thescattering wavefunction near the target, thereby diluting the impact of a resonance ina particular channel.
This chapter addresses the phenomenon of shape resonance in high harmonic
generation. It uses novel doubly adiabatic potential energy curves which illuminate theeffects of a strong field on energy barriers in particular channels[14]. These curves arethen used to construct a basis for the propagation of the time-dependent Schrödinger
equation. Finally, the effects of a shape resonance on a high harmonic spectrum areshown by comparing high harmonic spectra which are affected by a shape resonance tospectra which are not affected.
In the single active electron approximation, the Hamiltonian for an electron in thepresence of a spherically symmetric atomic potential U(r) and an external laser electricfield of strength F (t) is given byĤ = Tˆr + L̂2 /(2r 2 ) + V (r) + F (t)r cos(θ) = Tˆr + Ĥ|(r,t)(2.1)

Time dependent adiabatic channel functions and potential energy curves may befound by diagonalizing Ĥ|(r,t) at every value of (r, t)Ĥ|(r,t) φµ (r, t, Ω) = Uµ (r, t)φµ (r, t, Ω).

(2.2)

We perform this diagonalization in a spherical harmonic basis
φµ (r, t, Ω) =

lX
max
l=0

Aµl (r, t)Yl0 (Ω).

(2.3)

7
0

l=0
l=1
l=2
l=3
l=4
l=5
l=6
µ=0, E0=0.02 a.u.
µ=1, E0=0.02 a.u.
µ=2, E0=0.02 a.u.
µ=3, E0=0.02 a.u.
µ=4, E0=0.02 a.u.
µ=5, E0=0.02 a.u.
µ=6, E0=0.02 a.u.

V(r) (hartree)

-1

-2

-3

-4

0

5

10
r (bohr)

15

Figure 2.1: Adiabatic potential energy curves for Cs+ , calculated for E0 = 0 and E0 =.02. The potential curves depart from each other at long range, but are largely unaffectedin the region of the f-wave potential energy well.

Because cylindrical symmetry is not broken, m = 0 always.
In the limit that F (t) → 0, the problem reduces to the familiar case of an atom
in isolation, with Aµl (r, t) → δl,µ and Uµ (r) → ( l(l+1)+ V (r))δl,µ .
2r 2
The resulting potential energy curves serve to give insight into whether the shaperesonance survives the presence of an intense electric field. Figure 2.1 shows the potential
energy curves for Cs+ at various electric field strengths, and at various points in thelaser cycle.
As time goes by, this picture is complicated by the changing F (t), which has
the effect of mixing the different channels together. The time dependent Schrödinger
equation is written as an infinite, truncatable set of coupled differential equations[i∂/∂t − T̂r − Uµ ]Fµ (r, t) = −X1[ Qµν (r, t) + Wµν (r, t) + Pµν (r, t)∂/∂r]Fν (r, t) (2.4)
2v∂where Pµν =<< φµ (r, t, Ω)| ∂rφν (r, t, Ω) >> is the radial nonadiabatic coupling and∂Wµν =<< φµ (r, t, Ω)| ∂tφν (r, t, Ω) >> is the temporal nonadiabatic coupling.

By

the Hellman-Feynman theorem, Pµν = 3F (t)Cµν (r, t), and Wµν (r, t) = r Ḟ (t)Cµν (r, t),

8
where
Cµν (r, t) = <<φµ (r,t,Ω)| cos(θ)|φν (r,t,Ω)>>
Uν (r,t)−Uµ (r,t)µ 6= ν

0

µ=ν




(2.5)

The coupling due to the electric field is given by F (t)Zµν , where Zµν = r(Uν (r, t)−Uµ (r, t))Cµν (r, t). It is worth noting that the Cµν (r, t) term which appears in the electricfield and the nonadiabatic couplings will be finite for finite r. At zero field, Uµ (r) →00(l +1)( l(l+1)+ V (r))δl,µ and the energy denominator Uµ (r, t) − Uν (r, t) → δl,µ δl0 ,ν l(l+1)−l
.
2r 2
2r 2

When the field is nonzero, Uµ (r, t) and Uν (r, t) for µ 6= ν are distinct eigenvalues of theirreducible (since m is restricted to be 0) Hamiltonian. Thus Uµ (r, t) 6= Uν (r, t) and theresulting nonadiabatic couplings are finite.
In order to avoid time-dependent channel functions and couplings, time propagation was performed in the basis
χα (r, Ω) =XFµ(α) (r)φµ (r, Ω)

(2.6)

µ

of eigenfunctions of the Hamiltonian calculated with the field at maximum strength andsetting Wµν = 0. Here the time-independent Schrödinger equation is given by(Eα − T̂r − Uµ )Fµ(α) (r, t) = −X 1( Qµν (r) + Pµν (r)∂/∂r)Fν(α) (r, t),
2
ν

(2.7)

where Eα is the energy eigenvalue corresponding to basis state χα (r, Ω).
This calculation was performed using FEM-DVR basis functions[45]. These functions are orthonormal in the limit of Gauss-Lobatto quadrature and yield a bandedmatrix for the kinetic energy operator T̂r . Herman-Skillman atomic potentials [22] wereused for V (r). In addition, a complex absorbing potential was added in order to absorb
flux past a cutoff point. The resulting Hamiltonians were therefore sparse and complex
symmetric but not Hermitian. Use of the non-Hermitian Hamiltonian requires use ofa different orthogonality relation: hφ∗µ | |φν i = δµ,ν , where the extra conjugation on theleft-hand function causes the unconjugated function to appear in the integral[50].

9In the following calculations, an initial wavefunction was found by diagonalizing
the field-off Hamiltonian, then projected onto a basis of eigenfunctions of the field-on
Hamiltonian, calculated at the maximum of the electric field strength. For the fieldstrengths (E0 ≤ 0.02) used in these calculations, the field-on basis was truncated using
an energy cutoff of 4 hartree. As this energy is far above the harmonic and above
threshold ionization cutoffs for these fields, such a cutoff should not affect the calculated
spectra, and changing the cutoff from 4 to 8 hartree did not affect the results.
High harmonic spectra were calculated by taking the Fourier transform of thetime-dependent recombination dipole d(t) = hψg | z |ψ(t)i. The resulting HHG spectra
for hydrogen is compared with results found in [28] in Figure 2.2, which was found using
a different form of the dipole d(t) = hψ(t)| z |ψ(t)i. The disagreement is of the sameorder as was found in [28] to arise between the different forms of the dipole.
Several difficulties arise when trying to extract scattering information from an
observable such as the high harmonic spectrum. These difficulties arise from the complex complex rescattering wavefunction, which must undergo a complicated evolution
between the initial state and recollision. In order to contribute to d(ω), an electron
must tunnel free of the parent ion and propagate in the continuum before colliding once
again with the parent ion. Furthermore, there are two times during any laser cycle when
the original tunneling may occur, and these times will vary as the intensity envelope
increases or decreases. As a result, the recolliding wavefunction will have a complex
nature which will not be easy to separate from the high harmonic spectrum in order tofind the recombination dipole. An additional practical problem arises from the comblikenature of the spectrum: information from most scattering energies will be lost due todestructive interference between photons produced during different half-cycles of thelaser pulse. Thus, scattering properties reconstructed from the harmonic spectrum will
have data points separated by 2ω, making it difficult to resolve narrow scattering features. For a typical titanium-sapphire laser with driving wavelength 800 nm, ω = 0.0565

100.00012

2|d(ω)| (bohr )Current work
Krause et al

1e-08

1e-12

0

0.25

0.5
ω (Hartree)

0.75

1

Figure 2.2: Comparison between recombination dipole calculated in this work and by
Krause et al [28].

11
hartree, making the spacing between adjacent harmonics equal to 0.113 hartree = 3 eV.
For this reason, it would be desirable to perform such experiments using as small a valuefor ω as possible.
In order to minimize the effects of the continuum wavefunction, it is desirable tocompare two harmonic spectra with similar rescattering wavefunctions, but which result
from recombination into different orbitals. The effects of the rescattering wavefunction
can then be minimized by examining the ratio of the spectra as a function of therescattering electron energy Escat = w − Ip.
As will be shown in Chapter 4, the structure of the rescattering wavefunction is
primarily determined by the local potential felt by the electron at the time of rescattering; the impact of the initial orbital is mostly felt in the amplitude of the rescattering
wavefunction and an overall phase. For this reason, high harmonic spectra were calculated for two exited states of Cs+ : 5p5 6s1 and 5p5 5d1 , hereafter referred to as thes state and d state, respectively. The potential experienced by these electrons wasapproximated by a Herman Skillman potential which gave the bound state energy forthe d orbital as 0.24 H, compared with 0.30 for the s orbital. Due to the similarity inbinding energy, the ionization amplitude will be approximately the same as a functionof electric field strength for the two electronic states.
These two orbitals were chosen with the goal of seeing the f-wave shape resonanceknown to exist in xenon-like atoms [54, 55, 56]. Due to dipole selection rules1 , theresonant f-wave continuum electron can recombine into the d orbital, but not the s.
In the vicinity of xenon, both the d- and f-wave potentials rapidly become more
attractive with increasing nuclear charge. The f-wave potentials for a variety of atoms
in this region are shown in figure 2.3. By the time that Z=57 for lanthanum, the f-wave
potential has become sufficiently attractive to support a bound state. Cs + was chosen
1

Strictly speaking, l is no longer a good quantum number in the presence of the laser’s electric field.
However, near the atom the channel wavefunctions ψµ (r, Ω) are nearly diagonal in l, so arguments basedon selection rules are useful, if no longer exact.

12
4V(r) (hartree)

2

0
Cd
In
Sn
Te
I
Xe
Cs
Ba
La
Cs+

-2

-4

-6
0

1

2

r (bohr)3

4

5

Figure 2.3: In the vicinity of xenon, the f wave potential falls rapidly with increasing
nuclear charge. At lanthanum (Z=57), the potential supports a bound state. Here fwave Herman-Skillman potentials, are plotted for a variety of xenon-like atoms.

for the combination of having a relatively deeply bound 5d orbital at 0.23 H and a
low energy f-wave resonance, which occurs at a scattering electron energy of 0.1 H. In
addition, the relatively weak potential barrier may be affected more strongly by thelaser electric field than are the more robust barriers seen in Xe or Cs, thereby increasing
the visibility of electric field effects.
The photoionization cross section at zero field for the Cs+ 5d electron is shownin figure 2.4a, while figure 2.4b shows the p- and d-wave quantum defects. The shaperesonance is visible as a peak in the cross section, combined with a Π phase shift in the fwave quantum defect. Figure 2.5 shows the same information for the 6s photoionization
cross section, which does not exhibit a shape resonance.
As an additional aide to understanding the high harmonic spectra, photoionization calculations were made using field-on adiabatic potential energy curves, but neglecting nonadiabatic coupling terms in r and t. The associated angular functions wereapproximated as pure Yl,0 spherical harmonics. The resulting one-dimensional continuum wavefunctions were then normalized by matching the calculated wavefunctions to13

2

Photoionization Cross Section (bohr )

a)

0.1

0.05

0

0

Scattering Phase Shift

1.57

0.2

b)

0.8

1

0.8

1

f-wavep-wave

0

-1.57

0.4
0.6
Electron Energy (hartree)

0

0.2

0.4
0.6
Electron Energy (hartree)

Figure 2.4: a) Photoionization cross section vs energy for the 5d subshell of Cs + in the
5p5 5d1 excited state. b) p- and f-wave scattering phase shifts vs energy. The shaperesonance corresponds with a change in the f-wave phase shift of π.

142

Photoionization Cross Section (bohr )

0.005

0.004

0.003

0.002

0.001

0

0.2

0.4
0.6
Electron Energy (hartree)

0.8

1

0.4
0.6
Electron Energy (hartree)

0.8

1

b)

0

Scattering Phase Shift

a)

-0.1

-0.2

-0.3

-0.4

-0.5

0

0.2

Figure 2.5: a) Photoionization cross section vs energy for the 6s subshell of Cs + in the
5p5 6s1 excited state. b) p-wave scattering phase shifts vs energy for Cs+ . In contrast
to the 5d shell, no shape resonance effects are present.

15
to WKB solutions0

ψ(r ) =s2
(cos δ cos(πk(r0 )

Z r0

k(r)dr) + sin δ sin(

Z r0k(r)dr)).

(2.8)

Although nonadiabatic coupling is necessary for a complete calculation, the single
channel calculation illustrates how the relatively slight change to the potential energy
curves caused by a nonzero electric field can affect the cross section. Single channel
photoionization cross sections are shown for a variety of electric field strengths in figures
2.6 and 2.7. Figure 2.7 uses the peak electric field strengths experienced in the pulse,
while figure 2.6 uses lower electric field strengths, consistent with the strength of the
field at the time of recollision in the three step model, where the electrons return to theparent ion while the field is near a minimum. Figures 2.6c and 2.7c show the ratios ofthe d- and s cross sections as a function of the electron scattering energy, the quantitythat comparison of high harmonic spectra is intending to find.
Inspecting figures 2.6 and 2.7, it is apparent that the d state photoionization
cross section is relatively unaffected by the small electric field strengths used in thesecalculations. At higher field strengths, the maximum splits into three separate peaks.
In contrast, the s state photoionization cross section experiences a much greater change,
developing several maxima spaced by approximately 0.2 hartree. At the lower electricfield strengths consistent with the field strength at the time of recombination, almost
all variation in the d/s comparison spectra comes from changes in the s cross section,
while the d cross section undergoes relatively little change.
This pattern is repeated in the d/s ratios as interpolated from the high harmonic
spectra shown in figure 2.8. These spectra were made by first integrating the harmonic
intensity from (n − 1/2)ω to (n + 1/2)ω. The logarithms of the resulting intensities
were then interpolated as a function of the scattering electron energy E scat = nω − Ip tofind the envelope of the high harmonic spectrum as a function of energy. The ratio of
the recombination dipoles for the d and s spectra was then found by taking the square16

2

Photoionization Cross Section (bohr )

0.2

a)

0.15

0.1

0.05

0

0

0.2

2

Photoionization Cross Section (bohr )

0.02

0.4
0.6
Electron Energy (hartree)

b)

0.8

1

E0=.000 a.u.
E0=.001 a.u.
E0=.002 a.u.
E0=.003 a.u.
E0=.004 a.u.
E0=.005 a.u.
E0=.006 a.u.
E0=.007 a.u.
E0=.008 a.u.
E0=.009 a.u.
E0=.010 a.u.

0.015

0.01

0.005

0

40

d Cross Section / s Cross Section

E0=.000 a.u.
E0=.001 a.u.
E0=.002 a.u.
E0=.003 a.u.
E0=.004 a.u.
E0=.005 a.u.
E0=.006 a.u.
E0=.007 a.u.
E0=.008 a.u.
E0=.009 a.u.
E0=.010 a.u.

0

0.2

0.4
0.6
Electron Energy (hartree)

0.8

c)

1

E0=.000 a.u.
E0=.001 a.u.
E0=.002 a.u.
E0=.003 a.u.
E0=.004 a.u.
E0=.005 a.u.
E0=.006 a.u.
E0=.007 a.u.
E0=.008 a.u.
E0=.009 a.u.
E0=.010 a.u.

30

20

10

0

0

0.2

0.4
0.6
Electron Energy (hartree)

0.8

1

Figure 2.6: a) Cs+ 5d single channel photoionization cross section vs energy at low
electric field strengths. b)Cs+ 6s single channel photoionization cross section vs energy
at low electric field strengths. c) ratio of 5d to 6s photoionization cross sections vsenergy at low electric field strengths.

17a)E0=.000 a.u.
E0=.010 a.u.
E0=.011 a.u.
E0=.012 a.u.
E0=.013 a.u.
E0=.014 a.u.
E0=.015 a.u.
E0=.016 a.u.
E0=.017 a.u.
E0=.018 a.u.
E0=.019 a.u.

2

Photoionization cross section (bohr )

0.1

0.05

0

0

0.4
0.6
Electron energy (hartree)

0.8

0.05

0

d Cross Section / s Cross Section

30

0

0.2

0.4
0.6
Electron energy (hartree)

10

0

0.8

1

E0=.000 a.u.
E0=.010 a.u.
E0=.011 a.u.
E0=.012 a.u.
E0=.013 a.u.
E0=.014 a.u.
E0=.015 a.u.
E0=.016 a.u.
E0=.017 a.u.
E0=.018 a.u.
E0=.019 a.u.

c)

20

0

1

E0=.000 a.u.
E0=.010 a.u.
E0=.011 a.u.
E0=.012 a.u.
E0=.013 a.u.
E0=.014 a.u.
E0=.015 a.u.
E0=.016 a.u.
E0=.017 a.u.
E0=.018 a.u.
E0=.019 a.u.

b)

2

Photoionization Cross Section (bohr )

0.1

0.2

0.2

0.4
0.6
Electron energy (hartree)

0.8

1

Figure 2.7: a) Cs+ 5d single channel photoionization cross section vs energy at fieldstrengths comparable to pulse peak b)Cs+ 6s single channel photoionization cross section
vs energy at field strengths comparable to pulse peak strengths. c) ratio of 5d to 6s
photoionization cross sections vs energy at field strengths comparable to pulse peak.

18roots of the interpolated intensities as a function of scattering energy. Because theseratios are calculated from harmonic peaks spaced by 2ω, they cannot resolve narrow
features of the comparison spectra. As in the photoionization cross sections, there isa broad minimum at 0.1 H, followed by a maximum at 0.2H and associated with a πphase shift. At higher intensities, a second peak is seen at 0.5 H. This peak is harder
to interpret. At low field strengths, the d/s ratio as calculated in the single channel
photoionization develops a second peak at 0.6 H, slightly higher than the 0.5 H seenin the interpolated spectra. Such a disparity might arise due to nonadiabatic coupling,which is included in the time dependent calculation but not in the photoionization
cross section. Alternatively, the high harmonic peaks in this energy region, which are
beginning to approach the harmonic cutoff, may give a deceptive picture of the d/s ratio
in this energy region.
Overall, the d/s comparison spectra show a much stronger resemblance to photoionization calculations made at lower field strengths (figure 2.6c), than at higher fieldstrengths (figure 2.7c). This is consistent with the three step model, in which recollision
happens when the field is near a zero. The comparison spectra do show an observable dependence on electric field strength, however. This could indicate that electricfield does play an observable role in the recombination step of high harmonic generation.
Such an effect could also result if the flux of the rescattering wavefunctions at particularenergies had a different dependence on intensity for the s- as opposed to d- states.

2.1

Summary and Conclusions
This chapter has investigated the question of atomic shape resonances in an atom

subjected to an intense laser field. The overarching questions to be answered are whether
resonant features survive the introduction of the field, how these features are altered,and whether these features could be seen in an experiment.
Doubly adiabatic potential curves were found to give valuable insight for answer-

19

d cross section / s cross section

0.4

a)

E0=0.011 a.u.
E0=0.012 a.u.
E0=0.013 a.u.
E0=0.014 a.u.
E0=0.015 a.u.
E0=0.016 a.u.
E0=0.017 a.u.
E0=0.018 a.u.

0.3

0.2

0.1

0

0

0.2

6.282

0.4
0.6
Electron Energy (hartree)

0.8

1

0.4
0.6
Electron Energy (hartree)

0.8

1

b)

3.141

d phase - s phase0-3.141

E0=.011 a.u.
E0=.012 a.u.
E0=.013 a.u.
E0=.014 a.u.
E0=.015 a.u.
E0=.016 a.u.
E0=.017 a.u.
E0=.018 a.u.

-6.282-9.423-12.564

00.2

Figure 2.8: a) ratio of interpolated d harmonic spectrum to interpolated s harmoniccross section, normalized by ADK ionization rate b) difference in phase between d ands harmonic spectra.

20
ing these questions. Using these curves, it becomes easy to see whether features suchas a potential barrier survive the presence of the field. In addition, although the full
Hamiltonian includes nonadiabatic coupling between channels, single channel photoionization calculations were found to have a great deal of explanatory power for qualitative
features of the high harmonic spectrum.
The high harmonic spectra were found to have information about the scattering
wavefunction encoded in the harmonic peak heights, although this information is hidden by the complicated form of the scattering wavefunction. Much of this extraneous
information can be removed by taking the ratio of the spectra for multiple electronic
states of the same atom, since the rescattering wavefunction is largely independent of
the initial orbital from which the electron tunnels. These comparison spectra were
found to change with the applied field strength in a way which is qualitatively similar
to the single-channel photoionization calculations, although the correspondence is notexact. The comparison spectra were most similar to the photoionization calculations
when the electric field used in the photoionization calculations was considerably weaker
than the peak field strength of the pulse, consistent with a recombination occurring neara minimum of the electric field. The shape resonance itself was found to be relatively
insensitive to the field strength, although the potential barrier in Cs+ is not particularlyhigh. The peak of the resonance did not shift appreciably in photoionization calculations until E0 u 0.015. In a high harmonic experiment, where recollision usually occurs
near a minimum of the field, it is to be expected that shape resonance effects will still
be visible at much higher peak field strengths.

Chapter 3

Molecular scattering and the Limits of the Plane Wave
Approximation in the Measurement of Molecular Properties

As seen in the prior chapter, scattering phenomena such as shape resonances
can have measurable effects on scattering observables such as the harmonic spectra
measured in an experiment. This in turn has implications for the many current efforts tomeasure molecular properties by monitoring the scattering output channels in recollisionexperiments [32, 43, 41]. The best known such technique is molecular tomography[24], which uses high harmonic spectra from aligned molecules to reconstruct molecular
wavefunctions.
Rescattering electrons offer clear advantages as probes of molecular structure.
The intrinsic timescale of an ionization-acceleration-rescattering process is on the orderof a single half-cycle of the driving laser field, typically a few fs. Because the liberated
electron is accelerated by the driving laser field, simple formulas arising from classicalphysics are sufficient to map the energy of a rescattering event to the instants in a laser
half-cycle when the electron is liberated and returns, allowing time resolution to bepushed to the the sub-fs level.
Because such techniques all rely on the same underlying physical process wherein arescattering electron interacts with the parent ion, they are all inherently limited by the
degree to which the rescattering process is understood. However, to date most efforts
to measure molecular properties have treated the rescattering wavefunction as a free

22electron plane wave, unperturbed by the electron interaction with the parent ion. Priorwork relating to such reconstructions has dealt with bandwidth limitations arising fromthe HHG spectrum [39], orthogonality of the scattering- and bound-state wavefunctions[47, 42], and perturbative treatments of the ionic Coulomb potential [53]. This chapter
investigates the departure from plane wave scattering which is caused by a nonzero
molecular potential and the implications of that departure for molecular tomography.
“Tomographic images” of bound states are calculated for a one dimensional square well,and for two molecules in three dimensions, N2 and F2 . Much of the material presented
in this chapter appeared first in [66].

3.1

Scattering States and Ramifications for Molecular Tomography
At its heart, the tomographic procedure attempts to measure the dipole matrix

element
d~~k =

Zd3 ~xψ~k (~x)~xψg (~x)(3.1)

in momentum space between a continuum wavefunction ψ~k (~x) which asymptotically~goes as eik·~x and a particular orbital ψg (~x) of some target molecule. In the limit where
the molecular potential is zero, the plane wave approximation for the scattering states
would be exact, and the wavefunction could be reconstructed according to~xψg (~x) =Z~

d3~ke−ik·~x d~k(3.2)

A nonzero molecular potential complicates this picture. In one dimension, the
WKB approximation gives the continuum scattering state as

where k(x) =

Rx
1
0
0
ψc ∝ pei k(x )dx
k(x)

(3.3)

p
2(E − V (x)). In the vicinity of the molecule, both the amplitude andthe phase of the scattering state depart from the plane wave approximation.

23
An ideal tomographic experiment would measure d~~k between continuum states
and an unperturbed molecular ion. In contrast, in rescattering experiments, recombination occurs in the presence of a strong and time-varying external laser field. The
magnitude of the incoming wavefunction is affected by tunnel ionization from the highest
occupied molecular orbital (HOMO) and the propagation of the electron between ionization and recombination. In addition, the high harmonic spectrum is sharply peakedat frequencies which are multiples of the driving laser frequency. For these reasons,
it would be very appealing to measure d~~k using photoionization rather than high harmonic generation [47]. It is not clear how the phase of d~~k would be measured in such an
experiment, but it is at least conceivable to do so by introducing some kind of beamsplitter with interfering pathways. However, since this chapter is concerned with thelimitations to tomographic reconstruction under ideal circumstances, the photoionization amplitudes d~~k will be treated as though they could in principle be found. The
issues discussed here with respect to tomographic reconstruction for a photoionization
experiment apply also to HHG tomography, with the stipulation that the relevant scattering states should be calculated in the presence of an external laser field in order toprovide an exact description.
Problems with the tomographic reconstruction procedure arise when the scattering states ψ~k (~x) begin to deviate from the plane waves that were assumed in the initial
theoretical formulations [24]. As can be seen in equation 3.3, this deviation becomespronounced when the potential experienced by the electron is comparable to the scattering energy. In this case, the measured d~~k will depart from the Fourier transform of~xψg (~x).
In equation 3.1, substitution of
ψ~k (~x) = (2π)

−3/2

Z

~0
d3~keik ·~x ψ̃~k (k~0 )

(3.4)

24
and evaluating the integral over d3 ~x yields
d~~k = (2π)−3Z^
d3 k~0 ψ̃~k (k~0 )(~xψg )k~0(3.5)^
where (~xψg )(k~0 ) represents the Fourier transform of ~xψg (~x), the quantity which tomographic procedures hope to measure, and ψ̃~k (k~0 ) represents the Fourier transform of
the scattering state ψ~k (x). In Equation 3.5, the scattering states define a Fourier-space^
mapping from the desired function (~xψg )(k~0 ) to the measured function d~~k . In general,
this mapping will not be diagonal, as the scattering states ψ~k (~x) will have Fourier components at k~0 6= ~k due to distortions by the molecular potential. This mapping is not
generally invertible without knowledge of the molecular scattering states.
A 1D square well provides a simple example whose study can document theextent to which the electronic potential energy affects the outcome of a tomographic
reconstruction based on the plane wave approximation. A potential of the form
 V |x| ≤ x0
V (x) =
 0 |x| > x0

(3.6)

yields scattering states

where k2 =Aeikx + Be−ikx x ≤ −x0(scat)
ψ|k| (x) =Ceik2 x + De−ik2 x |x| ≤ x0






Eeikx + F e−ikx
x > x0(3.7)√
k 2 − 2V .

The two linearly independent solutions ψ±|k| (x) are now chosen such that theiroutgoing wave components go as e±ikx as x → ±∞. For ψ|k| , this corresponds toA=

4e2i(k+k2 )x0 kk2,−(k − k2 )2 + e4ik2 x0 (k + k2 )2
B=0C=

2ei(k+3k2 )x0 k(k + k2 )−(k − k2 )2 + e4ik2 x0 (k + k2 )2

(3.8)
(3.9)
(3.10)

25
D=−

2ei(k+k2 )x0 k(k − k2 )−(k − k2 )2 + e4ik2 x0 (k + k2 )2E=1(3.11)(3.12)

F =e2ikx0 (e4ik2 x0 − 1)(k 2 − k22 )−(k − k2 )2 + e4ik2 x0 (k + k2 )2(3.13)

A=

e2ikx0 (e4ik2 x0 − 1)(k 2 − k22 )−(k − k2 )2 + e4ik2 x0 (k + k2 )2(3.14)

and for ψ−|k| ,

B=1
C=−
D=

2ei(k+k2 )x0 k(k − k2 )−(k − k2 )2 + e4ik2 x0 (k + k2 )22ei(k+3k2 )x0 k(k + k2 )−(k − k2 )2 + e4ik2 x0 (k + k2 )2
E=0

F =

4e2i(k+k2 )x0 kk2−(k − k2 )2 + e4ik2 x0 (k + k2 )2(3.15)(3.16)(3.17)(3.18)(3.19)

Note that although k(x) takes on only two values in this problem, each scatteringsolution has nonzero Fourier components for k 0 6= k, k2 .
The true dipole matrix elements may now be compared to the plane wave approximation. A two-node bound state wavefunction was chosen as a simple example which
nonetheless possesses nontrivial spatial structure. Setting x0 = 2.5 bohr, V = −1.61
hartree yields a two-node wavefunction with E = −.5 hartree. Figure 3.1a compares
dipole matrix elements calculated between this bound state and continuum functions
described by scattering eigenfunctions and plane waves. Whereas the plane wave dipole
matrix elements are all purely imaginary, calculating the matrix elements using scattering states gives both real and imaginary components. If the dk calculated using
scattering states are now treated as “measured” dipoles and used to construct a tomographic image of the original bound state using equation 3.2, the resulting image will
be complex valued. Figure 3.1b compares the original bound state with its tomographic
image.

26

2a)

Re(scattering)Im(scattering)Re(plane wave)
Im(plane wave)

dk (bohr)

1

0

-1

-2
-10

0.6

-5

0
-1
k (bohr )

5

10

b)

0.4

ψ(x)

0.2
0
-0.2
-0.4
True bound state wavefunction
Re(tomographic image)Im(tomographic image)

-0.6
-10

-5

0
x (bohr)

510Figure 3.1: The tomographic reconstruction procedure applied to the 1D square well.
(Top) Comparison of dipole matrix elements dk = hψc (x)| x |ψg (x)i, calculated using
plane waves and scattering states for ψc (x). (Bottom) Because the scattering statematrix elements differ from those calculated using plane waves, the reconstructed imageof the orbital will differ from the true bound state wavefunction.

27
The principles seen in the case of the 1D square well also limit tomographic reconstruction in true molecular systems, although molecular systems are much more computationally challenging due to the complicated potentials which describe the electron-ion
interaction. We calculate the electron-ion scattering states using FERM3D [57], a code
which is designed for the highly non-centrosymmetric potentials seen in these systems.
This potential is described by Vmol = Vs + Vex + Vpol , where Vs is the local electrostatic potential, Vex is the exchange potential arising from antisymmetrization of the
wavefunction and treated in the local density approximation, and Vpol is a polarization
potential that describes the relaxation of the target under the influence of the incoming
electron. Figure 3.2 compares the total photoionization cross section for N 2 calculatedwith FERM3D and the plane wave approximation to a prior calculation and experiment[6], while Figure 3.3 compares cross sections for F2 to a previous calculation[27]. For
both molecules, FERM3D gives cross sections with sizes comparable to prior calculations, with photoionization maxima shifted higher than in the comparison. In both
molecules, the plane wave cross sections are too large by a factor of five.
The calculated dipole matrix elements may also be used to tomographically reconstruct the molecular orbitals. In a tomography experiment, the experimentally measurable quantity is the dipole matrix element between a bound state of the molecule~and a scattering state whose incoming-wave portion asymptotically goes to e ik·~x .
For incoming-wave boundary conditions, FERM3D calculates dipole matrix elements(−)

dqE;l,m = hψE;l,m | ²̂q · ~x |ψg i(3.20)(−)

where ψE;l,m is an energy-normalized wavefunction which obeys incoming-wave boundary conditions [57](−)

lim ψE;l,m =r→∞X

l0 ,m0Ylm (r̂)(2i)−1 (fl+ (r)δl,l0 δm,m0 −†fl− (r)Sl,m;l0 ,m0 ).

(3.21)

28
30

FERM3D
Plane Wave / 5
Carravetta
Experiment

Photoionization Cross Section (Mb)

25

20

15

10

5

0

20

30

40
Photon Energy (eV)50

60

Figure 3.2: N2 photoionization cross sections vs photon energy. Calculations madeusing FERM3D and the plane wave approximation are compared to experimental measurements taken from [6].

Photoionization Cross Section (Mb)

30

20

10

0

FERM3D
Plane Wave / 5
Kilcoyne, GIPM/D
Kilcoyne, GIPM/A
20

40
Photon Energy (eV)60

Figure 3.3: F2 photoionization cross sections vs photon energy. Calculations made using
FERM3D and the plane wave approximation are compared to theoretical calculations
taken from [27].

29
where f ± (r) has the asymptotic formf ± (r) → (πk)−1/2 r±i/k e±ikr±iη , ² = k 2 /2 ≥ 0(3.22)

are radially outgoing/incoming Coulomb spherical waves, as defined in [4] and σ l =arg[Γ(l + 1 − i/k)] is the Coulomb phase shift.
To find d~q , it is now necessary to find the superposition ψ~k = Alm (~k)ψE;lm whosek~outgoing component matches the outgoing component of eik·~x . Expanding~eik·~x = 4π

X∗il jl (kr)Yl,m (x̂)Ylm(k̂)(3.23)

l,m

where jl (kr) →r→∞ (2i)−1 (ei(kr−lπ/2) − e−i(kr−lπ/2) ) are spherical Bessel functions,matching coefficients of Ylm (r)eikr yields∗
Alm (~k) = 4πil e−iσl Ylm(k̂)k 1/2

(3.24)

where the factor of k 1/2 converts the energy-normalized matrix elements calculated in
FERM3D to momentum normalization.
d~~k is now given by
d~~q =k

X q
d~l,m (k 2 /2)Al,m (~k)(3.25)

l,m

As in the 1D square well, the tomographic image of the orbital may now be computed by substituting the calculated d~~k for the experimentally measured quantity in the
tomographic reconstruction procedure. As photoionization is not limited to the molecular HOMO, tomographic images may be calculated for all the orbitals of a molecule.
Such tomographic images will in general be complex-valued, and will differ according
to which polarization component is used for the tomographic procedure.
We present tomographic images of various orbitals of N2 and F2 , calculated in
the body-fixed frame using the x,y, and z polarization components. Each image is given
with the real and imaginary components, and includes an orange bar (not always visible)
extending between the two atoms of the molecule for the purposes of scale. For N 2 ,30
tomographic reconstructions for the 1πu (Fig. 3.4), 3σg (Fig. 3.5), 2σu (Fig. 3.6) and
2σg (Fig. 3.7) orbitals are shown. For F2 , reconstructions were calculated for the 1πg(Fig. 3.8), 3σg (Fig. 3.9), 1πu (Fig. 3.10), 2σu (Fig. 3.11), and 2σg (Fig. 3.12) orbitals.
For these example molecules, tomographic reconstruction tends to preserve theσ or π, gerade or ungerade character of the orbitals in question. However, the reconstructed orbitals may display additional radial nodes not found in the original orbitals.
Features which correspond to features of the original orbitals may be distorted in shapeand size, and display a spatially varying complex phase. Finally, tomographic imagesof the same orbital made using different polarization information may produce differing
images of the same orbital. Many of these features are also seen in [42], which treats thescattering process from a multielectron perspective, but does not consider the distorting
effects of a molecular potential.

3.2

Conclusions
The use of rescattering electrons as a probe of molecular properties offers many

exciting avenues for future research. However, the rescattering process is itself more
complicated than has been recognized in early reconstruction efforts, and is worthy of
study in its own right.
For molecular tomography, the results presented in this paper suggest that at
energy scales where such distortion is significant, tomographic reconstructions may besignificantly distorted from the “true” orbitals these methods seek to find. Nevertheless,for the example molecules presented here, the tomographic reconstruction procedure wasable to successfully reproduce the σ or π, gerade or ungerade nature of the orbitals inquestion. Reconstructions made using differently polarized dipole matrix elements gavedifferent tomographic images of the same orbital, while reconstructions made from aparticular polarization gave tomographic images with spatially varying complex phase.
Both of these properties could be useful as an experimental check of reconstructed31

Figure 3.4: Comparison of the Hartree-Fock orbital and associated tomographic images
for the N2 1πu orbital. a) and b) give the real and (zero) imaginary components of the
Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized dipole matrix element. e) and f) give the realand imaginary components of the tomographic image made from the y-polarized dipole
matrix element. g) and h) give the real and imaginary components of the tomographic
image made from the z-polarized dipole matrix element.

32

Figure 3.5: Comparison of the Hartree-Fock orbital and associated tomographic images
for the N2 3σg orbital. a) and b) give the real and (zero) imaginary components of the
Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized dipole matrix element. e) and f) give the realand imaginary components of the tomographic image made from the y-polarized dipole
matrix element. g) and h) give the real and imaginary components of the tomographic
image made from the z-polarized dipole matrix element.

33

Figure 3.6: Comparison of the Hartree-Fock orbital and associated tomographic images
for the N2 2σu orbital. a) and b) give the real and (zero) imaginary components of the
Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized dipole matrix element. e) and f) give the realand imaginary components of the tomographic image made from the y-polarized dipole
matrix element. g) and h) give the real and imaginary components of the tomographic
image made from the z-polarized dipole matrix element.

34

Figure 3.7: Comparison of the Hartree-Fock orbital and associated tomographic images
for the N2 2σg orbital. a) and b) give the real and (zero) imaginary components of the
Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized dipole matrix element. e) and f) give the realand imaginary components of the tomographic image made from the y-polarized dipole
matrix element. g) and h) give the real and imaginary components of the tomographic
image made from the z-polarized dipole matrix element.

35

Figure 3.8: Comparison of the Hartree-Fock orbital and associated tomographic images
for the F2 1πg orbital. a) and b) give the real and (zero) imaginary components of the
Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized dipole matrix element. e) and f) give the realand imaginary components of the tomographic image made from the y-polarized dipole
matrix element. g) and h) give the real and imaginary components of the tomographic
image made from the z-polarized dipole matrix element.

36

Figure 3.9: Comparison of the Hartree-Fock orbital and associated tomographic images
for the F2 3σg orbital. a) and b) give the real and (zero) imaginary components of the
Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized dipole matrix element. e) and f) give the realand imaginary components of the tomographic image made from the y-polarized dipole
matrix element. g) and h) give the real and imaginary components of the tomographic
image made from the z-polarized dipole matrix element.

37

Figure 3.10: Comparison of the Hartree-Fock orbital and associated tomographic images
for the F2 1πu orbital. a) and b) give the real and (zero) imaginary components of the
Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized dipole matrix element. e) and f) give the realand imaginary components of the tomographic image made from the y-polarized dipole
matrix element. g) and h) give the real and imaginary components of the tomographic
image made from the z-polarized dipole matrix element.

38

Figure 3.11: Comparison of the Hartree-Fock orbital and associated tomographic images
for the F2 2σu orbital. a) and b) give the real and (zero) imaginary components of the
Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized dipole matrix element. e) and f) give the realand imaginary components of the tomographic image made from the y-polarized dipole
matrix element. g) and h) give the real and imaginary components of the tomographic
image made from the z-polarized dipole matrix element.

39

Figure 3.12: Comparison of the Hartree-Fock orbital and associated tomographic images
for the F2 2σg orbital. a) and b) give the real and (zero) imaginary components of the
Hartree-Fock orbital. c) and d) give the real and imaginary components of the tomographic image made from the x-polarized dipole matrix element. e) and f) give the realand imaginary components of the tomographic image made from the y-polarized dipole
matrix element. g) and h) give the real and imaginary components of the tomographic
image made from the z-polarized dipole matrix element.

40
wavefunctions: ideal reconstructions would have a spatially uniform complex phaseand reconstructions made with different polarization information should agree with one
another. Additionally, experiments could use higher scattering energies to minimize thescattering state distortions due to interactions with the molecular potential.
The sensitivity of the scattering states to the molecular potential also offers theprospect for new types of experiments. Such experiments could monitor the movementof charge within the parent ion at ultrafast timescales. For example, a two-centerinterference experiment of the type discussed in [32] might observe the movement ofcharge in a diatomic molecule by observing interference maxima/minima occurring atdifferent energies for the short and long rescattering trajectories.

3.3

Appendix: Gauges and Dispersion Relations
Within the overall framework of the plane wave approximation, several heuristic

methods have been suggested to improve the accuracy of tomographic reconstructions[32]. For the 1D square well, we considered the effects of phenomenological “dispersion
relations” and reconstructions made in the momentum, rather than velocity gauge.
A dispersion relation attempts to correct for an electron’s shorter wavelength byp~
substituting ei~q·~x for eik·~x in Equation 3.2, where |~q| = 2(k 2 /2 − ²V ), where V is thepotential felt by the electron in the interaction region, and ² ∈ [0, 1].

Tomography may also be performed in gauges other than the length gauge given
in equation 3.2. If both continuum and bound wavefunctions are eigenstates of the
Hamiltonian, the dipole matrix element is identical in the length(l)dk = hψg | x |ψk i(3.26)

hψg | ip |ψk i
Eg − E k(3.27)and velocity
(p)

dk =

gauges. From the momentum gauge form of the dipole matrix element, and employing

41
the plane wave approximation for |ψk i, it is possible to generate a second tomographic
reconstruction
ψ

(p)

(x) =

Z

2
(p) Eg − k /2

dqeiq(k)x dk

−q(k)

.

(3.28)As a plane wave is not an eigenfunction of the scattering Hamiltonian, this reconstruction will in general give a different image of the target orbital than a reconstruction
made using the length gauge.
We tested both the length- and momentum-gauge tomographic reconstructions
pusing dispersion relations q(²) = 2(k 2 /2 − ²V ), using the overlap of the true wavefunction and its (normalized) tomographic image as a figure of merit. The test employed thesame V = −1.61, x0 = 2.5 potential and E = −.5 target wavefunction used to generate
Figure 3.1.
Figure 3.13a gives the magnitude of the overlap between the two tomographic
images and the ground state and between each other as a function of the dispersion
parameter ². For this choice of potential and target orbital, the tomographic reconstructions gave a very poor overlap with the target orbital for ² ∈ [0, 1], reaching amaximum magnitude of 0.60 at ² = 0.26 in the momentum gauge. In the dipole gauge,
the maximum overlap was achieved at ² = 1.0, also giving an overlap of magnitude 0.60.
Figure 3.13b compares the maximally overlapping reconstructions to the true groundstate wavefunction.
A perfect reconstruction would give the same image regardless of the gauge thetomographic procedure was performed in. However, agreement between images made inseparate gauges does not guarantee the accuracy of the reconstruction. Although bothgauges gave nearly identical tomographic images at ² = 0, the resulting images gaveamong the worst overlaps with the target orbital.

42

1

a)

Overlap of length gauge image with true wavefunction
Overlap of momentum gauge image with true wavefunction
Overlap of length gauge image with momentum gauge image

0.8

Overlap

0.6

0.4

0.2

0

0

0.2

0.4
0.6
Dispersion parameter ε

0.8

1

b)
0.4
0.2

ψ(x)

0
-0.2
-0.4

True bound state wavefunction
Re(length gauge image)Im(length gauge image)Re(momentum gauge image)Im(momentum gauge image)

-0.6
-0.8
-20

-10

0
x (bohr)

1020

Figure 3.13: a) Magnitude of the overlap between (normalized) tomographic images ofa bound wavefunction and the true
p wavefunction, and between different tomographic
images, calculated using q(²) = 2(k 2 /2 − ²V ), ² ∈ [0, 1]. b) Comparison of the maximally overlapping tomographic images to the true wavefunction. In the momentum
gauge, maximal overlap was obtained for ² = .26, while in the length gauge, maximal
overlap was obtained for ² = 1. Both images have been normalized and rotated to givea purely real overlap with the true wavefunction.

Chapter 4

High Harmonic Generation in SF6 : Raman-excited Vibrational
Quantum Beats

4.1

Introduction
In the prior chapters, and indeed in most treatments of high harmonic generation,

HHG is treated as a purely electronic process. A single active electron tunnels free ofa parent molecule and eventually scatters from it, but the molecule itself is treatedin essentially the same way as a lone atom, as no more than a complicated potential
influencing the active electron. In reality, molecules differ from atoms in an essential way
because they possess non electronic internal degrees of freedom, which can themselves
be affected by the high harmonic process.
The possibility that such internal degrees of freedom could play a detectable rolein HHG is intriguing, because the intrinsic timescale in HHG – the time necessary toionize, propagate and rescatter – is only half a laser cycle. This is faster than many
chemically interesting processes, holding out the possibility that HHG could serve asa probe of molecular motion. Alternatively, tailoring the state of a molecule prior to
HHG could serve to give additional control over the generated light.
These issues were brought to the fore by an experiment at JILA [63]. In theexperiment, a high harmonic generating laser pulse was preceded by a weaker pulse
whose effect was to stimulate Raman-active vibrations in SF6 molecules. Varying thedelay between the two pulses had the effect of modulating the intensity of the HHG light44
generated by the second pulse. Moreover, the modulation occurred at frequencies equal
to those of the Raman-active normal modes stimulated by the first pulse. Although SF 6
has 3 non-Raman-active vibrational modes, none of their frequencies were detected in
the modulated signal.
This chapter presents a fully quantum mechanical model of high harmonic generation in molecules. This model provides a framework to interpret the observed modulation of high harmonic intensities observed in the JILA experiment, and is easily
extended to systems with more complicated dynamics. Secondly, it presents a versionof the three step model which has been improved for the purpose of treating HHG in
molecules with relevant internal degrees of freedom. Finally, the modulations predictedby this improved model are compared with the modulations observed in the JILA experiment. This chapter recapitulates and extends work which originally appeared in
[65].

4.2

The Vibrational Wavefunction of the Molecule
An important difference between atomic and molecular systems is that molecules

have nontrivial structure in the form of vibrational degrees of freedom. For an M atom
molecule, this corresponds to N = 3M − 6 (N = 3M − 5 for linear molecules) internal
degrees of freedom, which can be expressed in normal mode coordinates. The vibrational
wavefunction of the molecule can then be expanded as the product of simple harmonic
oscillator basis functions in each of the normal modes of the molecule:|ψvib i =∞Xn1 ,n2 ,...nN =0

An1 ,n2 ,...,nN |n1 , n2 , . . . nN i

(4.1)

where
a(i)† a(i) |n1 . . . nN i = ni |n1 . . . nN i

(4.2)

where |n1 n2 . . . nN i is the outer product of simple harmonic oscillator state |n1 i in thefirst normal mode, |n2 i in the second normal mode, and so on.

45
For the purposes of this chapter, all operators will be expanded to first orderin raising a†(i) and lowering a(i) operators in each normal mode i. At this level ofapproximation, the evolution of the vibrational wavefunction becomes separable, and theoverall coefficient An1 n2 ...nN (t) can be factored into the product of individual coefficientsani (t) and an overall function of time f (t)An1 n2 ...nN (t) = an1 (t)an2 (t) . . . anN (t)f (t).

(4.3)(Note that the coefficients ani (t) should not be confused with lowering operators a(i) .)In the interests of simplicity and clarity, the remainder of this chapter will use a
1-D picture, describing the evolution of the vibrational wavefunction for a single normal
mode. The concepts from the 1-D model extend readily to the nonseparated problem
in the event that operators involving two or more raising or lowering operators become
important. 1

4.3

Vibrational Interference
Although high harmonic generation is primarily an electronic process, the vibra-

tional state of the molecule can affect the harmonic intensity. This occurs because at
several points in the high harmonic process, the vibrational wavefunction has amplitudes either to stay unchanged or to “hop” up or down from simple harmonic oscillator
state |ni to states |n ± 1i, much in the same way that a photon in a beamsplitter has
amplitudes to take two or more paths. As in a beamsplitter experiment, two or moreindistinguishable pathways interfere with one another and modulate the output signal
detectably. The multiple pathways at play in the high harmonic process are diagrammedin Figure 4.1.
Raman Excitation

The first opportunity to change vibrational states occurs

when the vibrationally cold molecules are subject to the weak, off-resonant initial pulse.
1

The lack of modulation at combination frequencies like 2ω1 , ω1 + ω2 ,ω1 − ω2 , etc. suggests but
does not prove that such higher-order terms are not important in discussing the JILA experiment.

46

|n=0>

|n=1>

a0

a1e-iwt
M
I
N

e-

R
d0

d1

Figure 4.1: The vibrational interference model [65] in one dimension. The molecule ends
the first (Raman) pulse in a superposition of the ν = 0 and ν = 1 vibrational states.
After a time delay, the two vibrational states are mixed by stimulated Raman scattering(transfer matrix M), “hopping” during ionization (I) and recombination (R), as well as
evolution of the ionic wavefunction while the electron is away (N). Interference between
adjacent vibrational states modulates the high harmonic signal.

47This causes the molecule to undergo stimulated Raman scattering. Beginning the first
pulse in the |0i vibrational state, it is driven into a coherent superposition of the zerothand first vibrational state|ψvib i = a0 (0) |0i + a1 (0) |1i(4.4)

by the end of the pulse. (For the pulse length and intensity used in the JILA experiment,
calculations show no appreciable population of the |2i or higher states after the Raman
pulse. Accordingly, the |2i and higher states have been dropped from this analysis.)
Because only one normal mode is used in this analysis, the N quantum number has
been dropped, so that a0 = a10 (0) and a1 = a11 (0). t = 0 is chosen at some time after
the end of the first pulse, when the stimulated Raman scattering is over.
The vibrational state coefficients follow equations of motion given byµ¶
11X
iȧni (t) = ωi ni +ani (t) −EA (t)EB (t)×
2
2
A,B

£

√

αAB ani + ∂i αAB ( ni + 1ani +1 +

√

(4.5)

¤

ni ani −1 ) .

Here ωi is the normal mode frequency, indices A and B run over {x, y, z}, EA (t) isthe component of the electric field in the (body-frame) A direction at time t, Q i isthe normalized displacement associated with normal mode i and αAB (Q1 , Q2 , ...) is
the polarizability tensor of the molecule. These equations of motion have off-diagonal
elements only if ∂i αAB ≡ (2mωi )−1/2 ∂αAB /∂Qi |Qi =0 6= 0, which is the condition for amode to be Raman active. The polarizability tensor and its derivatives are found by
performing an unrestricted Hartree-Fock calculation [16] using the aug-cc-pVTZ basis
set. [11]Between the two pulses, the |0i and |1i states evolve as eigenstates of the simple
harmonic oscillator Hamiltonian. The states are then mixed once again by stimulated
Raman scattering during the high harmonic generating pulse. These effects are approximated by a unitary 2 × 2 transfer matrix M , where Mij is the amplitude to be in state|ii at the instant of ionization after beginning the second pulse in state |ji.

48
Ionization and Recombination

The molecule’s vibrational wavefunction

evolves further during each of the three steps –ionization, propagation, and recombination–of the three step model. The molecule’s vibrational state hops up or down a level during
ionization, evolves while the propagating electron is away from the molecule, and hops
once again when the electron recombines with the parent ion.
These hopping amplitudes arise because ionization and recombination, commonly
thought of as purely electronic processes, are both strongly modulated by molecular distortions. This is the simplest way in which the internal degrees of freedom in molecules
allows for behavior that has no analogue in atomic systems. Nonzero derivatives of ionization and recombination amplitudes translate directly into amplitudes for the molecule
to change its vibrational state during these processes. Taylor-expanding the ionization
operator about the equilibrium configuration of the neutral molecule,ˆˆ eq + ∂ I Q + O(Q2 ),Iˆ = I|
(4.6)
∂Q
√
∂ Iˆ
using the identity Q = (a + a† )/ 2mω and substituting I0 = î|eq , I1 = (2mω)−1/2 ∂Q,
the ionization operator can be rewrittenIˆ = I0 + I1 (a + a† ).

(4.7)

Identical logic gives the recombination dipole vector operator~ˆ = R~0 + R~ 1 (a + a† )
R

(4.8)

Vibrational Dynamics of the Parent Ion

Between the times of ionization

and recombination, the evolution of the internal state of SF+
6 is quite complicated.
This is because SF6 has three degenerate orbitals at the point of maximum symmetry.
Thus, at any nuclear configuration near this maximum symmetry point, these threeSF+
6 orbitals are very nearly degenerate, and are mixed with one another strongly bymolecular distortions. Orbital degeneracies can be broken and orbital energies cancross even with relatively small distortions. Because of this, it is necessary to treat

49
this interplay between electronic and vibrational states when describing the dynamics
of SF+
6 in the vicinity of the maximum symmetry point.
In its maximum symmetry configuration, SF+
6 belongs to the Oh point group,with three degenerate T1g orbitals which transform like axial vectors x̂, ŷ and ẑ. Whenthe molecule is distorted away from the maximum symmetry point via either an Eg ora T2g distortion, the triple degeneracy breaks up into three nondegenerate electronic
orbitals. The fully symmetric A1g or “breathing” mode preserves the triple degeneracy.
The full vibronic (vibrational-electronic) Jahn-Teller coupling matrix for a triply
degenerate system is given by[15, 38, 5]


√
gζ
gη 
g1 − gθ + 3g²


√


HT = gξ gζ
g1 − gθ − 3g²




gη
gξ
g1 + 2gθ(4.9)

This matrix represents the coupling between the states with x̂, ŷ, and ẑ symmetry,caused by vibrational operators, so that gζ represents the off-diagonal coupling between
the states with x̂ and ŷ symmetry, gη the off-diagonal coupling between the states withx̂ and ẑ symmetry, and gξ represents the off-diagonal coupling between the states withŷ and ẑ symmetry.
Defining the Eg normal mode coordinates Qθ and Q² , which transform like 2z 2 −x2 − y 2 and x2 − y 2 respectively, and the T2g normal mode coordinates Qξ , Qη and Qζ
respectively as the coordinates that transform like yz,xz, and xy, the functions g i are50
given by [15]
g1 =

1
1
2
2
2
2
2
2 kE (Qθ + Q² ) + 2 kT (Qξ + Qη + Qζ )(4.10)gθ =√1 VE Qθ + NE (Q2 − Q2 ) + N1 (2Q2 − Q2 − Q2 )
g
²
η
θ
ζ
ξ
3

(4.11)√
3N1 (Q2ξ − Q2η )g² =√1 VE Q² + 2NE Qθ Q² +
g
3

gξ =

√
VT2g Qξ + NT Qη Qζ + N2 Qξ ( 3Q² − Qθ )(4.12)(4.13)gη =√VT2g Qη + NT Qζ Qξ + N2 Qη (− 3Q² − Qθ )(4.14)gζ =VT2g Qζ + NT Qξ Qη + 2N2 Qζ Qθ

(4.15)

These constants were found by performing a CASSCF state-averaged calculation
for the three lowest energy states of SF+
6 for various displacements of the molecule awayfrom the maximum symmetry configuration. These energies are shown in Figure 4.2.
The CASSCF calculations were made using a basis of Hartree-Fock orbitals calculated
for neutral SF6 . These three adiabatic energies were then fitted to the eigenvalues ofthe diabatic Jahn-Teller coupling matrix. This process yielded VT2g =.001209 H/bohr,
VEg =.1406 H/bohr, N1 =-.0362 H/bohr2 , KT2g =.7288 H/bohr2 , KEg =1.8486 H/bohr2 .
For the A1g mode, which does not enter into the vibronic Hamiltonian, an adiabatic
potential E = VA1g QA1g + 1/2KA1g Q2A1g , with VA1g =.0645 H/bohr, KA1g =2.98 H/bohr2
gives the potential energy surface for all three electronic states.
Here there is a significant distinction between distortions of type Eg , which break
the triple degeneracy but have no linear off-diagonal terms, and distortions of typeT2g , which do contribute to off-diagonal coupling. For small distortions, an adiabatic
electronic state of an Eg distorted molecule will have the same symmetry – x̂, ŷ, ẑ – as the
diabatic electronic states. The adiabatic electronic states of a T 2g -distorted molecule,on the other hand, are linear combinations of the diabatic orbitals. An important
simplification is that VT2g , controlling off-diagonal coupling between different electronic
states, is small in SF+
6 and can be neglected for the short time between ionization andrecombination.

-993.33951a)E=E0+.064486 q+1.4911 q

2

Energy (hartree)-993.340-993.341-993.342-993.343-993.344-0.04-993.338

-0.02

b)

0
A1g displacement (bohr)

E=E0+0.14061 q + 0.8318 q

Energy (hartree)

0.04

0.02

0.04

0.05

0.1

2
2

E=E0- 0.14061 q + 0.8318 q

E=E0- 7.3653e-13 q + 0.11707 q-993.3400.02

2-993.342-993.344-993.346-0.04-993.338

-0.02

c)

0
Eg distortion (bohr)

2

E=E0+ 1.8342e-07 q + 0.21944 q
E=E0+ 0.0012091 q + 0.4369 q

2

2Energy (hartree)

E=E0- 0.0012091 q + 0.4369 q-993.340-993.342

-0.1

-0.05

0
T2g displacement (bohr)

Figure 4.2: a)Spherically symmetric A1g (breathing mode) distortions change electronic
state energies, but preserve the triple degeneracy. Nontotally symmetric E g (b) and T2g(c) distortions break the triple degeneracy of SF+
6 at the maximum symmetry point.
Adiabatic energies are fit to the eigenvalues of the vibronic coupling matrix (equation
4.9) to solve for the vibronic coupling constants.

52
The evolution of the ionic wavefunction is calculated in the vibrational basis of
the neutral molecule, but using potential energy surfaces calculated for the ion. Potential energy curves are found to quadratic order in Q using quantum chemistry calculations, then expressed in terms of raising and lowering operators by substituting
Q = (2mω)−1/2 (a + a† ), Q2 = (2mω)−1 (a + a† )(a + a† ). All terms up to linear in raising and lowering operators are then used to integrate the time-dependent Schrödinger
equation to find a transfer matrix N describing the evolution of the ionic wavefunction
between ionization and recombination.
Modulation of Harmonic Intensity

In the two-state model used here, thei−th vibrational wavefunction of the neutral molecule after recombination has occurredis |ψvib i = d0 |0i + d1 |1i, whereµd~0 d~1¶=µa0 (0)

a1 (0)e−iωτ

¶

T

~ .
M T IT N T R(4.16)

Here, e.g. M T denotes the transpose of matrix M .
The number of photons emitted in a given harmonic is proportional to d~0 · d~∗0 +d~1 · d~∗1 . The high harmonic intensity is a sum over all Raman active modes i:(i)

P (τ ) = P0 + Σi P1 cos (ωi τ + δi ) .

(4.17)The static P0 primarily results from terms of the form a0 (0)∗ a0 (0), while P1 results from~† ·terms of the form a0 (0)a1 (0)∗ eiωτ and a1 (0)a0 (0)∗ e−iωτ . Defining W = M † I † N † R~ IM , P0 = a0 (0)∗ W00 a0 (0) and P1 cos (ωt + δ) = 1 (a∗ eiωτ W10 a0 (0) + c.c.). Since I1
RN
2 1
and R1 are small relative to I0 and R0 , only their first-order terms are kept.

4.4

Describing the Continuum Electron
As seen in previous chapters, describing the evolution of the wavefunction for theactive electron is itself a challenging problem. When the electron first tunnels free ofthe parent molecule, its wavefunction is determined by both the molecular potential and53
the electric field of the laser. Once free, it propagates in the time-varying field of thelaser while feeling a weak force due to Coulomb attraction to the parent ion. Finally
it recollides with the parent ion, and is once again strongly distorted by the molecularpotential.
A full solution of the time dependent Schrödinger equation for this process would
be computationally demanding for complicated molecules such as SF6 . In addition,much of the information in the continuum wavefunction is not relevant to the HHG
problem: only a small part of the wavefunction overlaps with the unoccupied orbital
into which the rescattering electron recombines.
One frequently used treatment which avoids the complications of the full timedependent Schrödinger equation [7, 34] is based on a classical or semiclassical propagation of the continuum electron, ignoring the ionic Coulomb potential. The returning
electron wavefunction is then approximated as a plane wave throughout the recombination process. This approach has been successful in describing the high harmonic
cutoff, the chirp of the emitted high harmonic light and other quantities of interestin atomic systems [33, 26]. However, as was shown in an earlier chapter, the plane
wave approximation is not adequate to describe the returning electron because of thetremendous distortion caused by the electron’s interaction with the ionic potential andby exchange effects with the other electrons in the molecule. For the time-reversed problem of photoionization, it is known that the plane wave approximation is prone to error
for photoelectron energies smaller than the deepest k-shell binding energy. Energies
attained in high harmonic generation experiments usually fall below this range.
This section gives an improved semiclassical model of the free electron propagation, which can be connected to short-range calculations of the tunneling and rescattering wavefunctions using stationary phase techniques. Any method could be used
to calculate these short-range wavefunctions. In this treatment, the tunneling wavefunction is modeled semiclassically using ideas based on the initial value representa-

54
tion [37, 40]. The recolliding electron continuum wavefunction is described in terms ofelectron-molecule scattering states as used in chapter 3, calculated in the absence of an
external electric field.
In this way, the molecular potential affects the electron wavefunction at the twotimes when the electron is near the molecule. When the electron is far from the molecule,the comparatively simple evolution of its wavefunction is described using the shortesttime, dominant contribution to the Gutzwiller propagator [21]. Finally, stationary phasearguments serve to identify the isolated trajectories that encapsulate the effect of theelectron propagation in the field on high harmonic generation, greatly reducing thecomputational burden of propagating the continuum electron wavefunction.

4.4.1

Tunneling Ionization

During the ionization step, the tunneling electron wavefunction is described ina simple 1-D WKB tunneling picture, in which electrons are allowed to tunnel only indirections parallel to the laser’s applied electric field. This approach is motivated by
the semiclassical “Initial Value Representation” [37, 40], where a source wavefunction
acquires an imaginary phase (and hence an exponential rise or fall) along a trajectorythat passes through a classically forbidden region. The unperturbed molecular HOMO
here serves as the source wavefunction, so that the tunneling wavefunction is approximated by the unperturbed HOMO in the classically allowed region near the molecule,connecting to a WKB exponential which is set equal to the HOMO at the inner turning
point and decays exponentially until it reaches the outer turning point. SF6 has threedegenerate HOMOs: one of these is illustrated in Figure 4.3.
The classically forbidden region, illustrated in Figure 4.4, is defined by two turning
points, a and b, between which V (x) − E > 0. The slope of V (x) is C1 at inner turning
point a and C2 at outer turning point b. The tunneling wavefunction ψt (~r, t) is now
found by applying WKB connection formulas. The appendix derives the ratio of the55

Figure 4.3: One of three degenerate orbitals of SF6 . Red denotes positive lobes; blue
denotes negative lobes.

56

Figure 4.4: The tunneling wavefunction is approximated as an unperturbed molecular
HOMO inside the classically allowed region, connecting to a declining WKB exponential
in the classically forbidden region. Stationary phase trajectories leave from the outerturning point, beginning with zero velocity at time of ionization.

tunneling wavefunction at the outer turning point to the tunneling wavefunction at the
inner turning point
1
C1
Bi(0)
ψ(x = b)= e−Γ | |1/6
ψ(x = a)
2
C2
Ai(0)

(4.18)In this approximation, the tunneling wavefunction behaves like an Airy Bi function near the outer turning point, i.e. it has no linear complex phase term. This
property will be revisited in the section dealing with stationary phase analysis.

4.4.2

Semiclassical Propagation

When the active electron has tunneled free of the molecule, the evolution of its
wavefunction is controlled by the oscillating electric field of the laser, plus a residual
Coulomb attraction to the molecular ion. This relatively simple evolution continues until
the electron returns to the molecule, when the complicated molecular potential again
becomes significant. During this excursion, until it re-enters the non-Coulomb partof the potential, the continuum wavefunction can be approximated using Gutzwiller’s

57
semiclassical propagator[21]:X

K(~r, t; ~r0 , t0 ) =(2πi)−3/2

Classical TrajectoriespC(~r, t; ~r0 , t0 )×(4.19)

exp(iS(~r, t; ~r0 , t0 ) − iφ)R
Here S(~r, t; ~r0 , t0 ) is the action integral S = L(q, q̇, t)dt calculated for a classical2

Strajectory starting at (~r0 , t0 ) and ending at (~r, t), while C(~r, t; ~r0 , t0 ) = | −∂∂r0,A ∂rB |, whererA is the A-component of the vector ~r. φ is a phase factor equal to π2 times the numberof conjugate points crossed by the trajectory[21].
This propagator acts on the tunneling wavefunction ψt to give a semiclassicalcontinuum wavefunctionψc (~r, t) =

Z

d3~r0

Z

dt0 K(~r, t; ~r0 , t0 )ψt (~r0 , t0 ).

(4.20)

until the molecular potential asserts itself during the rescattering.
During the terminal portion of the scattering process, when the scattering wavefunction has a large dipole matrix element with the molecular HOMO to which it recombines, the electronic wavefunction is expanded into a truncated but in principle
complete basis of field-free electron-molecule scattering orbitals, calculated using techniques described in references [58, 59, 57]. Beyond the range of the molecular potential,
i.e. for r > r0 , the (l, m)-th independent scattering state is expressed as a partial wave±
expansion in terms of incoming and outgoing Coulomb radial functions fEl(r) and thescattering S-matrix as
1 −ψE,lm (~r) = √ fEl(r)Ylm (θ, φ)−
i 2
1 X +
√
f 0 (r)Yl0 m0 (θ, φ)Sl0 m0 ;lm (E), r ≥ r0 .
i 2 l0 m0 El(4.21)The laser electric field is typically far smaller when the electron returns to the ion. Neglecting its effect on the electron during its recollision with the ion, the time-dependent
wavefunction becomesψscat (~r, t) =

Z

dE

X
lm

Alm (E)ψE,lm (~r)e−iEt(4.22)

58
and the expansion coefficients Al,m (E) are given by
Alm (E) = e

iEt

Z

∗
d3~rψE,lm(~r)ψc (~r, t)(4.23)

for some chosen time t when ψc is projected onto the scattering states.
The dipole recombination amplitudes
dl,m (E) = hψE,lm | ²̂ · ~r |ψHOMO i(4.24)are calculated between the distorted scattering states and the molecular HOMO towhich the electron recombines.
Equations [4.20] and [4.23] now define the expansion coefficients Alm (E) for somechosen projection time t, in terms of two three-dimensional integrals over initial and final
positions, one integral over the time of ionization, and a summation over all possibleclassical trajectories. Stationary phase arguments dramatically simplify the calculationof these expansion coefficients.

4.4.3

Stationary phase calculation of scattering coefficients

The Gutzwiller propagator K(~r, t; ~r0 , t0 ) defined in equation [4.19] involves a summation over the classical paths having (~r0 , t0 ) and (~r, t) as their endpoints. Invoking theusual assumption that these paths are isolated, such classical paths have the property ofRextremizing the action integral S(~r, t; ~r0 , t0 ) = L(q, q̇, t)dt. If the initial or final pointsof the paths are varied, a new extreme-action path connects the new endpoints, and a
new extremal action corresponds to the new trajectory. From Goldstein et al [20] section
8.6, the change in the action integrals between the old classical path connecting the old

endpoints and the new path connecting the new endpoints is called the ∆ variation of
the action. The ∆ variation of the action is given by∆S(~r, t; ~r0 , t0 ) = (pA δqA − Ĥ(t)δt)|final
initial(4.25)

59
Because iS(~r, t; ~r0 , t0 ) appears in the exponent of the Gutzwiller propagator, slightchanges of ~r and ~r0 cause the phases of the integrals in equations 4.20 and 4.23 tovary rapidly, resulting in canceling contributions to the expansion coefficients. Thiscancellation occurs unless the ∆ variation of the action is counteracted by the changeof phase of the tunneling or scattering wavefunctions. A trajectory where the ∆ variation of the action is counterbalanced by the phase of the tunneling wavefunction at(~r0 , t0 ) and by the phase of the scattering wavefunction at (~r, t) is a “stationary phasetrajectory.” The initial and final points (~r0 , t0 ) and (~r, t) identify the points in the 7D
integral where the integrand oscillates slowly, giving a non-canceling contribution to the
expansion coefficients Alm (E).
To find such trajectories, expand the phase-oscillating parts of the 7D integral
Z
Z
Z
p
3
3
Alm (E) = d ~r d ~r0 dt0 eiEt C(~r, t; ~r0 , t0 )(4.26)∗
exp(iS(~r, t; ~r0 , t0 ) − iφ)ψE,lm(~r, t)ψt (~r0 , t0 ).

about starting point (~r0c , t0c ) and ending point (~rc , t)
The ∆ variation gives the expansion of the action integral
S(~rc + δ~r, t; ~r0c + δ~r0 , t0c + δt0 ) = S(~rc , t; ~r0c , t0c )+1 ∂2S
δrA δrB + p0A δr0A +2 ∂rA ∂rB
1 ∂2S
1 ∂2R
δr0A δr0B + Ĥδt0 +(δt0 )2 .
2 ∂r0A ∂r0B
2 ∂t20pA δrA +(4.27)

Finally, ignoring the angular derivatives of the spherical harmonics, the asymptotic phase evolution of the scattering states is given through the asymptotic phase±evolution of the Coulomb wave functions fEl
(r)
∗
ψE,lm
(~rc + δ~r) =(2πi)−3/2 −∗∗√ f (rc )Ylm
(θc , φc ) exp(ik(rc )δr)−−i 2 El

1
+
∗
√ fEl
0 Yl0 m0 (θc , φc )Sl0 m0 ;lm (E) exp(−ik(rc )δr)
−i
2
0
0
l ,m
X

where kl (r) =(4.28)p(2(E − Vl (r))). The tunneling wavefunction looks like a declining WKB

exponential in the forbidden region; it has no oscillatory phase part.

60
The condition for a nonoscillating integrand is now that the first-order terms inδ~r0 ,δ~r and δt0 must disappear. Setting the coefficients of δt0 to zero gives
p~20+ V (~r0 , t0 ) = EHOMO .
2m(4.29)

Setting the coefficients of δ~r0 to zero in Eq. 4.23 gives
p~0 = 0(4.30)and setting the coefficients of δ~r to zero in Eq. 4.23 gives
pr = ∓k(r)(4.31)

for the radial integral, where ∓ corresponds to the incoming and outgoing wave partsof the scattering states. Restricting consideration to trajectories that are incoming at
the time of projection, this last condition becomes
pr = −k(r) = −p2(E − V (r)).

(4.32)

Neglect of the angular derivatives of the spherical harmonics gives that the trajectorymust return with zero angular momentum.
Thus, a stationary phase trajectory is launched with zero momentum from theclassical turning point and returns to the molecule with zero angular momentum, andkinetic energy equal to the energy of the scattering state. This is a familiar result from,eg, [7, 33], with the distinction that the present work considers the effect of the electronion Coulomb interaction during the continuum propagation of the electron. Also, thewavefunction is projected onto scattering states shortly before recollision, rather thantreating the electron in a Volkov approximation throughout the recollision with themolecular ion. Because the electron-molecule scattering states are calculated with noexternal electric field present, there is a slight dependence on the time at which the
wavefunction is projected onto scattering states; here the projection is made when ωt

61for the laser cycle is equal to 3.9, i.e. when the short trajectories with energies equal tothe energy of the 39th harmonic have nearly returned to the molecule.
Because the linear phase variation of the integrand vanishes in the vicinity ofa stationary phase trajectory, the expansion coefficients Alm (E) are now found via
Gaussian integrals:i −∗∗Alm (E) = (2πi)−3/2 √ fEl(rc )Ylm(θc , φc )×2
exp(iS(~rc , t; ~r0c , t0c ) − iφ)ψt (~r0c )e−iEHOMO t0c ×Zi ∂2S(δt0 )2 ))×((I) = d(δt0 ) exp(
2 ∂t20Z∂2Si
δr0a δr0B ))×((II) = d3 (δ~r0 ) exp(
2 ∂r0A ∂r0B
Z
i ∂2S((III) = d3 (δ~r) exp(
δrA δrB ))2 ∂rA ∂rB

(4.33)

where the integrals labeled (I),(II) and (III) are evaluated as√∂2S
2πi| 2 |−1/2
∂t0
2
∂ S
|−1/2
(II) = (2πi)3/2 |∂r0A ∂r0B∂ 2 S −1/2|(III) = (2πi)3/2 |∂rA ∂rB(4.34)(I) =(4.35)(4.36)

yielding expansion coefficientsAlm (E) = (2πi)2
(

s

|

∂rA
|×
∂p0B

∂ 2 S −1/2 i
√ ×
)
∂t20
2

(4.37)−∗∗
fEl(rc )Ylm(θc , φc ) exp(iS(~rc , t; ~r0c , t0c ) − iφ)×ψt (~r0c )e−iEHOMO t0c
Once these expansion coefficients have been calculated, the dipole matrix element between the distorted scattering wave and the molecular HOMO is simply
~
D(E)
=

X
lm

Alm (E)d~lm (E).

(4.38)

62

4.5

Comparison with Experiment
The model of vibrational interference connects with this treatment of the high~ˆ is set to D(E).
~harmonic process when IˆRThis is broken up by setting Iˆ = ψt (~r0 , t0 ) and~ˆ = D(E)/~ˆ Both of these quantities are calculated for a molecule at the equilibriumRI.
geometry, and for a molecule displaced by 0.1 bohr in the normal mode coordinate. This
involves recalculating the scattering states and recombination dipoles for each distorted
molecular geometry.
The modulation of the 39th harmonic was chosen for purposes of comparison withexperiment, since this harmonic was considered in detail in [63]. The 39th harmonic
falls close to the measured cutoff, and can only be produced by a half-cycle coming close
to the maximum of the gaussian envelope of the laser pulse.
The JILA experiment used a gas jet as a source of SF6 , giving no preferred molecular orientation. Both ionization and recombination amplitudes are highly dependent onorientation. A rotational average was calculated for both the static and oscillatory partsof the harmonic intensity. Only those polarizations perpendicular to the propagating
laser beam for a given molecular orientation were included in these averages.
The non-totally symmetric Eg and T2g modes of the molecule are stimulated bydifferent amounts for different orientations of the molecule relative to the laser field. Ionization, recombination, and their derivatives are also highly orientation dependent. Because of this, although the Eg and T2g modes modulate the harmonic intensity strongly
for particular orientations of the molecule, the phase offset δ in P1 cos(ωτ +δ) from equation 4.17 changes with orientation, canceling some of the observed oscillation when the
angular average is performed. Figure 4.5 shows the spherically averaged peak-to-peak
modulation calculated when δ is allowed to vary for different molecular orientations andwhen δ is artificially set to 0 for all orientations. These two averages are comparedto the two experimental runs for which the modulations for all three modes could be63
distinguished from the background.
Although the agreement with experiment is not perfect, it is nevertheless significant that the simple model of vibrational interference presented here agrees withexperiment to the correct order of magnitude. This is particularly notable in light of
conventional Raman spectroscopy, in which the A1g peak is 20 times more prominent
than the others. It is difficult to precisely gauge the agreement of theory and experiment, due to the paucity of experimental data. Peak-to-peak modulations may varyby several percent from one experimental run to another [1] with SF6 , and at the timeof writing, only two experimental runs could distinguish all three Raman frequencies
from experimental noise. The modulation at 525 cm−1 , corresponding to the T2g mode,
appears most prominently in the experimental data, yet gives the smallest modulationin this treatment. The prominence of the T2g mode modulation may suggest that theoff-diagonal Jahn-Teller coupling VT2g is larger than given by these calculations. Alternatively, it may be necessary to model the experiment in more detail – spatially varying
laser intensity, uncontrolled carrier envelope phase, multiple laser half cycles, etc – thanincluded in the present work.

4.6

Conclusions
The problem of high harmonic generation in molecules may be conceptually sep-arated into two parts: the evolution of the continuum electron, and the evolution of theinternal (vibrational) wavefunction of the parent ion. This paper describes the evolution of the continuum electron in a model which combines a semiclassical treatment of
the propagation with a fully quantum mechanical description of the electron-molecule
scattering. This flexible and robust model has a simple conceptual link to existing semiclassical models, yet allows for a sophisticated treatment of the complicated electronmolecule scattering. The internal dynamics of the parent ion are tracked throughout thehigh harmonic process. Together, these two innovations serve to give an unprecedented

64

Experiment 1
Experiment 2
Phase information included
Phase information excluded-1

T2g mode, 525 cm-1

Peak-to-Peak Modulation

A1g mode, 775 cm
0.1-1

Eg mode, 642 cm
0.05

0
500

550

600

700
650
-1
Wavenumbers (cm )

750

800

Figure 4.5: Peak-to-peak modulation of the high harmonic signal vs. wavenumber,comparing theory to the two experimental runs for which data is available. For the bars
labeled “phase information included”, the angular average was performed allowing δ inequation 4.17 to vary as a function of angle. For the bars labeled “phase information
excluded,” δ was set to zero for all angles. Modulations corresponding to the same
frequency have been placed side-to-side for purpose of comparison. (Figure taken from
[65])

65
view of high harmonic generation in a large, complicated molecule with many internal
degrees of freedom, giving results which agree in order of magnitude to experiment.
The possibility that high harmonic generation may serve as an ultrafast interferometric probe of a molecular wavefunction is extremely promising. Such a wavefunctionneed not be prepared by an initial Raman pulse, as was done in the JILA experiment.
Instead, a preparatory pulse could photoionize a molecule, excite it to a higher electronic
state, or trigger the beginning of some other chemical process. In this way, wavepacket
evolution during chemical processes could be observed as it happens.

4.7

Appendix: Tunneling Ionization
In this approximation, the wavefunction in the forbidden region is found using

the WKB connection formulas. Illustrated in Figure 4.4, the classically forbidden region
is defined by two turning points, a and b, between which V (x) − E > 0. The slope ofV (x) is C1 at inner turning point a and C2 at outer turning point b.
Near turning point a, k 2 (x) ≡ 2m(E−V (x)) ≈ C1 (a−x) and the time independent
Schrödinger equation isψ 00 + k 2 (x)ψ = 0,(4.39)

which has solutions near x = a of
ψ(x) = Ai(

C1 (x − a)2/3
C1

)b1 + Bi(C1 (x − a)
2/3

C1

)b2

(4.40)

where Ai and Bi are Airy functions, asymptotically behaving likeAi(z) ⇒z→∞

3/2 ](2π)−1/2 z −1/4 exp[ −23 z(4.41)

Ai(z) ⇒z→−∞

π −1/2 (−z)−1/4 sin( 23 (−z)3/2 + π/4)(4.42)

Bi(z) ⇒z→∞

π −1/2 z −1/4 exp( 32 z 3/2 )(4.43)

Bi(z) ⇒z→−∞ π −1/2 (−z)−1/4 cos( 23 (−z)3/2 + π/4)
1/3

where z = C1 (x − a)(4.44)

66
Similarly, near x = b
1/3

1/3

ψ(x) = Ai(C2 (b − x))d1 + Bi(C2 (b − x))d2(4.45)

Under the barrier but away from the turning points, the WKB wavefunction is
given by
ψ(x) = π −1/2 |k(x)|−1/2 exp(

Z x1|k(x)|−1/2 exp(−2a|k(x0 )|dx0 ) sin(φ)−

Z x
a

0

(4.46)

0|k(x )|dx ) cos(φ)

for some value of φ.
Rb
Setting a |k(x0 |dx0 ≡ Γ, note thatZ b
Z x|k(x0 |dx0|k(x0 )|dx0 = Γ −(4.47)xaand for x close to b
Z x
b

2 1/2
k(x0 )dx0 = C2 (x − b)3/2
3(4.48)

Connecting the asymptotic forms of the Airy functions to the WKB solution in
the forbidden region gives the solution for x ≈ b as
−1/6

ψ(x) =x≈b 2C2

1/3

sin(φ)eΓ Ai(C2 (b − x))−1 −1/6
1/3
Ccos(φ)e−Γ Bi(C2 (b − x)),
2 2
giving a scattering phaseshift δ of
δ = tan−1 (4e2τ tan(φ))(4.49)(4.50)

yielding a resonance at φ ≈ 0, δ ≈ pi
2
Similar logic gives the wavefunction for x ≈ a−1/6

ψ(x) =x≈a C1

1/3

sin(φ)Bi(C1 (x − a))−(4.51)

1/6
1/3
C1 cos(φ)Ai(C1 (x − a))Finally, at resonance, φ = 0 and the ratio of the tunneling wavefunction at theouter turning point to the tunneling wavefunction at the inner turning point is
1
C1
Bi(0)
ψ(x = b)= e−Γ (| |)1/6
ψ(x = a)
2
C2
Ai(0)

(4.52)

Chapter 5

Laser-cluster interactions in the VUV energy range

Much of this thesis has dealt with the case of an isolated particle, exposed to alaser field which is treated classically. A single electron is liberated by the field, evolves
in the continuum, and interacts once again with the parent ion. In many ways, theinteraction of an atomic cluster with a vacuum ultraviolet (VUV) laser pulse does not
fit this scenario.
Intense lasers have only recently become available at these photon energies, andthe work done in this chapter was motivated by experiments [60] done at the free electron
laser (FEL) at the Tesla Test Facility (TTF) in Hamburg, Germany. The capabilities ofthis facility are documented in Ref. [17]. In this work, van der Waals clusters of xenon
atoms were exposed to 12.7 eV VUV photons, an energy range which had not previously
been explored.
What they found at this new energy range was very surprising. At a laser intensityof 7 × 1013 W/cm2 , isolated xenon atoms were found to produce only singly charged
ions [61, 48], corresponding to a single photon absorbed. However, when large clusters
were exposed to pulses of this strength, each atom was found to absorb up to 400 eV,corresponding to 30 photons! In addition, charge states of up to 8+ were detected.
At first glance, the reason for surprise might not be apparent. Previous experiments with noble gas clusters of krypton,argon or xenon atoms performed with infraredlasers (photon energies 1 eV) [23] with peak intensities on the order of 10 16 W/cm2

68
and pulse lengths on the order of 100 fs had shown very efficient absorption of laser energy. The atoms in the cluster ionize to form nanoplasmas, accompanied by high charge
states, strong x-ray emission and keV electrons[9, 51]. Experimental studies performed
by Lezius et al [35] with Ar and Xe clusters found that the nanoplasmas undergo very
energetic fragmentation, with ions of kinetic energies as large as 1 MeV ejected fromthe clusters.
In most models of the laser cluster interaction, designed to deal with infrared
driving lasers, the dominant heating mechanism is collisional heating[8, 10]. In sucha model, electrons tunnel free of their parent ions and are accelerated by the laser’s
electric field. However, because clusters have a very high local density of atoms, manyelectrons will scatter from atomic centers. It is therefore possible for an electron to beaccelerated one way by the laser’s electric field and scatter backwards from an atomic
center at just the right time to be accelerated in the opposite direction by the laser’s
electric field, which has switched directions in the meantime.
A back of the envelope calculation reveals that this can be a very efficient mechanism for heating in the case that the ponderomotive potential of the driving laser isvery large. For an applied electric field E = E0 cos(ωt), the kinetic energy of an electron
in one dimension may be given as
1
1
KE = (p − A)2 = (p2 + A2 − 2p · A)
2
2

(5.1)

where the vector potential A = − Eω0 sin(ωt) and the expectation value of p is < p >=E0ω sin(ωt).

If each backscattering changes the momentum from p to −p, the average
E2energy gain per collision is then 2 < p > A = 2 ω20 sin2 (ωt), which in turn averages to
E02= 4Up 1 .
ω2

In the case where the cluster is subject to intense infrared lasers, with strong
1

More calculations must be done to deal with the (much more likely) case of small angle elastic
scattering. However, as the change in energy comes from the p · A term, the answer will be proportional
to Up , just as in the example.

69
electric fields and small laser frequencies, Up can be very large and such a collisional
heating mechanism can deliver great amounts of energy to the cluster plasma. However,VUV photons have very high frequency, while available sources cannot produce laser
intensities as high as those which can be achieved with infrared lasers. In the Hamburg
experiment, the free electron laser was able to produce intensities of 7 × 10 13 W/cm2 ,for a ponderomotive potential of only 62 meV, nearly three orders of magnitude smaller
than the energy absorbed per atom!
Previous work [49] had indicated that inverse bremsstrahlung could serve as aneffective heating mechanism. In this picture, electrons freed from their parent ions by aninitial photoionization would form a small cluster plasma. They would be able to scatterfrom a xenon ion and absorb a photon from the laser field in the process (called inverse
bremsstrahlung scattering). Santra et al made the further approximation that inverse
bremsstrahlung cross sections should be calculated using realistic Herman-Skillman potentials rather than the pointlike Coulomb potentials used in prior treatments. Usingsuch a potential, an electron passing near an atomic center would feel an effective charge
Zef f much greater than the total charge of the ion, since the inner-shell electrons would
screen the nucleus less effectively. Because inverse bremsstrahlung cross sections go
as Z 2 , this could potentially have very large impacts on the rate at which the cluster
plasma could absorb photons.
In [49], Santra and Greene had found that using Herman-Skillman atomic potentials increased the number of photons absorbed per atom from about 1 (using Coulomb
potentials) to nearly 30. However, that study calculated inverse bremsstrahlung rates
using perturbation theory for both the electron-photon and the electron-ion interaction,and did not consider such effects as collisional ionization or recombination, or the finite
size of the clusters.
This chapter presents a model of the laser-cluster interaction in the VUV regime
which includes the effects of non-Coulombic potentials as well as the evolution of the70cluster as it absorbs energy from the laser. The electron-ion interaction is treated
nonperturbatively using variational R-matrix methods. The electron-photon interaction is treated using first-order perturbation theory. Although no electron correlation
is used, this treatment is expected to give more realistic photoionization and inverse
bremsstrahlung cross sections than those used in Santra et al. The model treats photoionization, collisional ionization and recombination, inverse bremsstrahlung heating,evaporation of electrons from the cluster, and expansion of the cluster due to hydrodynamic pressure of hot electrons and Coulomb repulsion. This work appeared first in
[64].

5.1

Photoionization
One distinction between the Hamburg experiment and previous experiments donewith IR lasers is that a single 12.7 eV photon has sufficient energy to ionize the highest
bound electron of xenon, which has an ionization potential of only 12.1 eV. Thus,
although the oscillating electric field is too weak to field ionize the xenon atoms, there
is still an efficient mechanism to create an initial cluster plasma.
The cross section for the transition from the bound state |φi i to the continuum
state |φf i is given by [19]σf i (E) = 4π 2 αω|~π · ~rf i |2 ,(5.2)

where α is the fine-structure constant, ω is the photon energy, ~π is the polarization
vector for the radiation, and ~rf i = hφf | ~r |φi i is the dipole matrix element coupling theinitial and final states of the electron. E is the kinetic energy of the photoelectron,
whose wave function in Eq. (5.2) is energy-normalized.
For linearly polarized light, chosen to be polarized in the ẑ direction, the relevant
matrix element is
hφf | z |φi i = IR (li , lf )ZdΩYl∗f mf (Ω) cos θYli mi (Ω) .

(5.3)

71
Here,
IR (li , lf ) =Z ∞drUf (r)rUi (r) ,

(5.4)

0

where U (r) = rR(r) denotes the rescaled radial wave function.
At a photon energy of 12.7 eV, only the 5p electrons of xenon can be ionizedby the laser field. Thus, the photoionization cross section may be calculated using
only electrons from this shell. Using an independent electron model to deal with the qelectrons in this shell, the total atomic photoionization cross section can be found by
averaging over the initial and summing over the final one electron states to give4αωσPI = q π 2
3 2li + 1
© 2ª
2
× li IR(li , li − 1) + (li + 1)IR(li , li + 1)

(5.5)

where li = 1 is the angular quantum number of the p shell. In calculating this result,the identities [46, 12]ZdΩYl∗1 m1 (Ω) cos (θ)Yl2 m2 (Ω) =



and



p
 l1
(2l1 + 1)(2l2 + 1)(−1)−m1 
−m1



1 l2  l1 1 l2 

 ,
0 m2
0 0 0











l1 + 1 if l2 = l1 + 1






X  l1l2
l3   l1
l2
l30  δ(l3 , l0 )δ(m3 , m0 )
3
3
,


=
2l
+
1
3
0
m1 ,m2 m1 m2 m3
m1 m2 m3

2

l 1 1 l 2 (2l1 + 1)(2l2 + 1) 
 = l1



0 0 0




0


if l2 = l1 − 1

(5.6)

(5.7)

(5.8)

otherwise

l2
l3 
 l1
were used.  represents the Wigner 3-j symbol, related to the Clebschm1 m2 m3Gordan coefficient by


l2
l
(−1)l1 −l2 −m l1√< l1 m1 l2 m2 |l1 l2 l − m > .
 =2l + 1m1 m2 m72
The radial integrals IR (li , lf ) [Eq. (5.4)] were calculated in the acceleration representation for wave functions generated by a variational eigenchannel R-matrix calculation [4] using a Herman-Skillman atomic potential [22].
Because of the efficient photoionization, the atoms in the cluster ionize quickly
to form a dense nanoscale plasma at an early stage in the laser pulse. This plasma has
the effect of screening the atomic potential felt by both bound and continuum electrons.
This lowers the ionization potential and changes both the initial- and final-state electron
wave functions. Because of this, cross sections for photoionization become larger as thescreening length in the plasma become shorter. With sufficient screening, it becomes
possible for ions to undergo additional photoionization.
To account for this process, the screened radial matrix elements were calculatedusing the same R-matrix methods as for isolated Xe atoms. However, before the initialand final-state wave functions were calculated, the Herman-Skillman potential was multiplied by a Debye screening factor exp (−r/λD ). (The electron Debye length is definedp
as λD = T /(4πne ) [29], where the electron temperature T in this expression is given inunits of energy. ne is the electron density.) The resulting matrix elements and ionization
potentials were then spline-interpolated in order to calculate the photoionization crosssection for a particular screening length. For most calculations, this screening length
was restricted to be no less than 4.64 bohr, the Wigner-Seitz radius of xenon at liquiddensity. Shorter screening lengths are discussed in a later section.

5.2

Inverse Bremsstrahlung Heating
A second effect of having a high density of free electrons in the cluster plasma isthat these electrons can undergo both stimulated and inverse bremsstrahlung processes,
creating a second mechanism by which energy can be absorbed in the cluster. Suchprocesses are the quantum mechanical analogue of the collisional heating mechanism
discussed earlier. Stimulated bremsstrahlung refers to photoemission into the laser

73
mode by an electron colliding with an ion in the plasma, while inverse bremsstrahlung
refers to the opposite process of an electron absorbing a photon in such a collision.
Because each bremsstrahlung event results in 12.7 eV deposited into or removed fromthe plasma, such processes serve a large role in heating the plasma.
Collisions of electrons with cluster ions are treated as independent events, allowing
heating rates to be found by calculating cross sections for a single electron-ion collision.
The cross section per unit energy for a free-free transition from initial state |φ E 0 ,l0 ,m0 i
to final state |φE,l,m i can be shown using Fermi’s golden rule, to equal¯¯2¯
4π 2 α ¯¯∂V0 ,l0 ,m0 i¯
σE,l,m←E 0 ,l0 ,m0 =
|φ
hφ
|
E,l,m
E
¯ .
ω3 ¯
∂z

(5.9)

E = E 0 − ω in the case of photon emission, E 0 + ω in the case of photon absorption.
Equation (5.9) describes the interaction of linearly polarized radiation in the acceleration
representation. V is the plasma-screened atomic potential experienced by the scatteredelectron.
As with photoionization cross sections, radial wave functions were calculated using a nonperturbative eigenchannel R-matrix approach. Matrix elements between theenergy-normalized wave functions were calculated in the acceleration gauge, where thewhere the 1/r 2 long-range dependence of ∂V /∂z ensures that the radial integral will
converge, although the continuum electron wave functions are not spatially normalizable.
Bremsstrahlung processes act as a powerful heating process even though microscopic reversibility ensures that cross sections for absorption and emission cross sections
have the same magnitude. This occurs because lower energy states are more highly populated than higher energy states in a thermal distribution. This study assumed thatthe electron gas reequilibrates rapidly after each bremsstrahlung event due to frequent
electron-electron collisions. Heating rates are then calculated assuming that the electron

74
probability distribution ρ(E) is at all times given by a Maxwell-Boltzmann distribution
ρ(E) = 2r

E −E/Te.
πT 3(5.10)

Using this distribution, it is possible to derive an expression for the change inelectron temperature with respect to time due to inverse bremsstrahlung processes.
The cross section defined in Eq. (5.9) describes a free-free transition between
orbital angular momentum eigenstates. We therefore introduce ρ(E, l, m), which is theprobability per unit energy to find an electron in the state |φE,l,m i. Clearly, ρ(E) =Pl,m ρ(E, l, m). If the wave function is normalized within a large sphere of radius R(not to be confused with the cluster radius), then the largest l that contributes to this√sum at a given kinetic energy E is lmax = R 2E [2]. Since, in thermal equilibrium,ρ(E, l, m) can depend only on E, we see that
ρ(E, l, m) =

ρ(E)
2R2 E

(5.11)

in the limit of large R (lmax À 1).
We are interested in radiation-induced heating, i.e. in the change of the electron
temperature due to photon absorption and emission. To this end, we will derive from
∂T
2
=
∂t
3

Z ∞
0

dEE

X ∂ρ(E, l, m)
l,m∂t(5.12)

a rate equation for the electron temperature, expressed in terms of the cross sections
for stimulated and inverse bremsstrahlung (Eq. (5.9)). As a final step in writing down
the equation for the time evolution of ρ(E, l, m), we must take into consideration that
ρ(E, l, m) refers to (spherical) box normalization, while the cross section per unit energy
in Eq. (5.9) is based on energy-normalized wave functions. For the sake of consistency,
it is necessary to change the initial state in the free-free radiative transition in Eq. (5.9)
from energy normalization to box normalization. This has the effect of multiplying the√cross section by π 2E/R.

75
Hence, if there are Na atomic scatterers within the normalization volume and alaser beam of intensity I, the rate of change of ρ(E, l, m) is given by√I 2π∂ρ(E, l, m)
= Na
×
∂t
ω R
Xn
√
√
σE,l,m←E−ω,l0 ,m0 E − ωρ(E − ω, l0 , m0 ) + σE,l,m←E+ω,l0 ,m0 E + ωρ(E + ω, l0 , m0 )l0 ,m0o√√− σE+ω,l0 ,m0 ←E,l,m Eρ(E, l, m) − σE−ω,l0 ,m0 ←E,l,m Eρ(E, l, m) .

The first row in the curly brackets in Eq. (5.13) describes the population of |φ E,l,m i via
photoabsorption (photoemission) from states with energy E − ω (E + ω); the second
row describes the depopulation of |φE,l,m i due to photoabsorption and photoemission
from this state. Equation (5.13) implies a nondegenerate electron gas.
An electron state with energy E will then communicate with states of energy
E − ω, which are on average more densely populated than itself. Since the absorptionand emission cross sections are equal, this tends to populate the state of energy E
and depopulate the states of energy E − ω, resulting in a net heating. The state will
also communicate with states of energy E + ω, which are less densely populated thanitself, thereby again tending to populate the higher-energy states while depopulating
the lower-energy state.
Combining Eqs. (5.10), (5.11), (5.12), and (5.13), we are led in a natural way tothe following definition of the inverse bremsstrahlung cross section (per unit energy):
σE+ω←E =XX

σE+ω,l,m←E,l0 ,m0 .

(5.13)

l,m l0 ,m0Using Eqs. (5.6), (5.7), (5.8), and (5.9), this can be written as
σE+ω←E =

4 2α
π
3 ω3
X©
ª
lJR2 (l, l − 1) + (l + 1)JR2 (l, l + 1) ,×

(5.14)lwhere
JR (l, l0 ) =Z ∞0drUE+ω,l (r)dVUE,l0 (r) .
dr(5.15)

76
The rate of change of the electron temperature due to inverse bremsstrahlung is then
∂T
∂t

=

µ ¶3/2 hi
22π
1 − e−ω/T
na I
9
T
Z ∞
×
dEe−E/T σE+ω←E .

(5.16)

0

The parameter na stands for the number of atoms per unit volume. In general, theions in the plasma are not all in the same charge state. Denoting the fraction of Xe i+
P
(i)
(i)
by f (i) , σE+ω←E in Eq. (5.16) is replaced with i f (i) σE+ω←E , where σE+ω←E is theinverse bremsstrahlung cross section in the field of Xei+ .

Figures 5.1 and 5.2 illustrate the dramatic effects of the ionic potential on the
inverse bremsstrahlung cross section. In Fig. 5.1, as the scattering electron collides
with the ion at higher and higher initial energies, it probes regions of the ionic potential
at which the ion nucleus is screened increasingly poorly by inner-shell electrons. As aresult, the inverse bremsstrahlung cross section rises to many hundreds of times that of
the naked Coulomb potential.
Adding plasma screening to this picture has the effect of supplementing the screening effects of inner-shell electrons with the screening effects of plasma electrons. As aresult, the scattering electron feels the effects of the ionic nucleus more strongly thanin the pure Coulomb case, but less strongly than in the case of the unscreened ionic
potential. This is seen in a steady decrease of the inverse bremsstrahlung cross section
as the screening range decreases.

5.3

Collisional Ionization and Recombination
Although photoionization and inverse bremsstrahlung are the only processes bywhich the cluster can absorb photons from the laser beam, they are not by themselves
enough to explain the cluster’s evolution. As the pulse progresses, large numbers of
free electrons fill the cluster. These electrons can liberate other electrons via collisional
ionization if they have sufficient energy, or they can undergo three-body recombination

Inverse Bremsstrahlung Cross Section (Mb/hartree)

77
Coulomb Potential
Herman-Skillman Potential
10000

1000

100

10

1

0

50
Initial Electron Energy (hartree)

100

Figure 5.1: Inverse bremsstrahlung cross sections [Eqs. 5.13 and 5.14] for an electronwith incident energy E to absorb a 12.7-eV photon are given for an electron in the fieldof a purely Coulombic 1+ potential and for an electron in the field of a Xe HermanSkillman atomic potential. The effects of atomic structure on inverse bremsstrahlung
rates are quite pronounced. Figure taken from [64].

Inverse Bremsstrahlung Cross Section (Mb/hartree)

12000

10000

8000

6000

4000

2000

0

0

50
Initial Electron Energy (hartree)

100

Figure 5.2: The inverse bremsstrahlung cross section as a function of energy is shown
for an electron in the field of a Debye-screened Xe Herman-Skillman potential, with
the Debye screening length λD ranging from 1 a.u. to 20 a.u. As λD grows, the crosssection approaches the limit of no plasma screening, shown in this graph by the dottedline. As the Debye length of the cluster plasma shrinks, the charged ion is shielded more
effectively from the scattering electron, and the inverse bremsstrahlung cross section isdecreased. Figure taken from [64].

78

15

Energy (hartree) per Atoma)Total Energy Absorbed
Electron Kinetic Energy
Energy Absorbed due to Photoionization

10

5

0

-200

-100
0
100
Time (fs) Relative to Pulse Maximum

200

1
b)

0+
1+
2+
3+

Fraction of Population

0.8

0.6

0.4

0.2

0

-200

-100
0
100
Time (fs) Relative to Pulse Maximum

200

Figure 5.3: Evolution of a 1500 atom cluster exposed to a 100 fs, 7 × 10 13 W/cm2
pulse, employing only photoionization and inverse bremsstrahlung heating. a) Energy
absorbed vs. time. b) Ionic population vs. time. Xe2+ and Xe3+ are produced efficiently
via photoionization. Figure taken from [64].

79with an ion, thereby heating the plasma.
These processes do not directly involve laser photons, but they nevertheless play
an important role in the heating of the cluster. Because inverse bremsstrahlung goes as
the square of the ionic charge, even small populations of highly charged ions can havegreat impacts on the rate at which the cluster is heated. In addition, the ability for acold plasma to gain energy through recombination and a hot plasma to lose it throughionization allows the cluster ions to act as an energy reservoir for the plasma electrons
throughout the pulse2 .
Including the effects of ionization and recombination, the rate equation for thenumber per unit volume ni of charge species i is given by
∂ni
∂t

=

I i−1
ni−1 − σPi I ni )(σω PI+Si−1 ni−1 ne − Si ni ne(5.17)+Ri+1 ni+1 n2e − Ri ni n2e ,where ne is the number of electrons per unit volume. The photoionization cross sections
i were calculated in Section 5.1. The ionization and recombination coefficients S andσPIiRi for the reaction Xei+ +e− →Xe(i+1)+ +2e− are calculated later in this section.
The evolution of the cluster is found by integrating these equations numerically,
along with equations for the energy in the free electron gas and the radius of the cluster.
As a general rule, this set of equations is quite stiff. We performed this integration using
the Rosenbrock method [44].
There are two requirements for a satisfactory treatment of collisional ionization
and recombination in the cluster. First, both processes must occur at appropriate
rates. Second, the rates for ionization and recombination must be consistent with one
another, in the sense of driving the cluster toward chemical equilibrium at all times. The
2

During the pulse, energy is typically deposited into the cluster too quickly for the charge state
distribution to equilibrate. Thus, although the charge state distribution may act as an energy reservoir
for the plasma electrons, it is typically far from being in thermodynamic equilibrium.

80
second requirement is particularly significant because the usual treatment of collisional
ionization (also used in this study) uses the semiempirical Lotz formula [36] for ionizationfrom the j th subshell,SijÃ
Z ∞ −x
j
a
q
e
1
j
−7
i
= 6.7 × 10
dx
3/2
Pj /T Pj /T x
T!
Z
bj exp (cj ) ∞
e−y
cm3
−
dy
,
Pj /T + cj Pj /T +cj ys(5.18)

to find ionization coefficients. In the Lotz formula, aj , bj , cj are semiempirical constants,
qij the number of equivalent electrons Xei+ contains in the j th subshell, Pj the ionization
potential in eV, and T the temperature in eV. For charge states of 0, . . . , 5+, semiempirical constants corresponding to the 5p sublevel are chosen. For charge states of 6+
and 7+, which have no 5p electrons in the ground state, constants are chosen whichcorrespond to the 5s sublevel. Because this formula is only a semiempirical approximation, it is important to use recombination coefficients which are consistent with theionization coefficients to prevent the model from settling into an incorrect equilibrium
charge distribution.
The ratio between ionization and recombination coefficients can be obtained using
the concept of equilibrium constants. In a plasma at equilibrium, the rate of collisions
ionizing Xei+ to form Xe(i+1)+ + e− must be equal to the rate at which Xe(i+1)+ + e−
recombines to form Xei+ . However, the relative populations of reactants and byproducts
at equilibrium is a thermodynamic property of the system which does not depend on
intimate knowledge of ionization and recombination rates. Recombination coefficients
consistent with the correct relative populations at equilibrium can then be applied tomodeling the cluster plasma, which is not in general in a state of chemical equilibrium.
These two rates are given respectively by
Si ni ne = rate of ionizing collisions(5.19)

81
andRi+1 ni+1 n2e = rate of recombining collisions .

(5.20)

Hence, the ratio S/R is given by
neq neqSie= i+1eq .
Ri+1
ni(5.21)eq
eq
The fraction (neqi+1 ne )/ni is known as the equilibrium constant for the reaction, andcan be calculated thermodynamically.
In any reaction A → B+C at equilibrium, the chemical potentials for the forwards
and backwards reactions must be balanced µA = µB + µC . The chemical potential ofeach species is given by a partial derivative of the Helmholtz free energy,
µi =

∂F
.
∂Ni T,V(5.22)The Helmholtz free energy is given by F = −T ln Ztot , where Ztot is the partition
function for the system as a whole.
Factoring the total partition function into the product of individual particle partition functions (which implies independent particles),Ztot (NA , NB , NC , V, T ) =ZA (NA , V, T )ZB (NB , V, T )ZC (NC , V, T )zA (V, T )NA zB (V, T )NB zC (V, T )NC
,
=
NA !
NB !
NC !

(5.23)

and assuming Ni >> 1 yieldsµi = −T∂ ln(Zi (V, T ))zi
= −T ln( ) .
∂Ni
Ni(5.24)

Imposing balanced chemical potentials yieldszB (V, T )zC (V, T )
NB NC
=
NA
zA (V, T )(5.25)

or equivalently
Keq (T ) =

nB nC
=
nA

zB (V,T ) zC (V,T )V
V,zA (V,T )V(5.26)

82
where the equilibrium constant Keq is a function of temperature only.
If the ionization potential of Xei+ is given by Pi , then the partition functions are
given by

zi+1
ze

Z ∞

dEe−E/T ρi (E) ,0
Z ∞−Pi /T= e
dEe−E/T ρi+1 (E) ,0Z ∞=dEe−E/T ρe (E) .

zi =(5.27)

0Through most of the lifetime of the pulse, tight plasma screening destroys the Rydberg
states and most of the internal degrees of freedom of the various ions, leaving the densityof states ρ(E) dominated by the center of mass term and by a combinatorial term m 
(5.28)
D(i) =  naccounting for the number of ways n electrons can be distributed in m orbitals. For

charge states up to 6+, we use m = 6, n = 6−i. For 7+ and 8+, we use m = 2,n = 8−i.
If we exploit this by setting ρi (E)/D(i) = ρi+1 (E)/D(i + 1), the common integral in zi
and zi+1 falls out of the equilibrium constant, yielding
3

(Xe
Keq

i+ →Xe(i+1)+ +e− )T 2 D(i + 1)
= e−Pi /T √ 3
2π 2 D(i)(5.29)

Equation (5.29) can now be combined with Eqs. (5.19) and (5.21) to yield recombination
rate coefficients which have appropriate magnitude and which, in combination with theionization coefficients, drive the system toward the correct equilibrium distribution at
all times.
A gas of charged particles has different thermodynamic properties from an idealgas due to Coulomb interactions between the constituent particles. Zel’dovich andRaizer [67] calculate the adjustment to Keq due to a Debye-Hückel potential. The
equilibrium constant including Coulomb effects can be written
32

T D(i + 1)
(Xei+ →Xe(i+1)+ +e− )Keq= e−(Pi +∆Pi )/T √ 3,
2π 2 D(i)(5.30)

83
where the change in ionization potential due to Coulomb effects is ∆Pi = −(Qi +1)/λD ,the Coulomb potential between the ion core and an electron held at distance λ D .
In our approach, by explicitly calculating bound state energies for Debye-screened
Hartree-Slater potentials, we calculate this adjustment to the ionization potential directly. Our adjustment behaves similarly to the Zel’dovich and Raizer correction, but
is larger for longer screening lengths and smaller at shorter screening lengths.
One advantage to the equilibrium constant approach is that it conceptually separates information about thermodynamic balance from the rate at which the systemseeks that balance. As a result, any formula for ionization or recombination coefficients
could be substituted for the Lotz formula, with the accuracies of the overall rate and of
the equilibrium constant used the only criteria for validity of the formula.
Including the effects of collisional ionization and recombination has a pronounced
effect on the evolution of the cluster. In Fig. 5.3, the evolution of the cluster is calculated
employing only photoionization and inverse bremsstrahlung. In contrast, Fig. 5.7 showsthe evolution of the same cluster employing photoionization, inverse bremsstrahlung,collisional ionization and recombination, and evaporation of energetic electrons fromthe cluster. Allowing ionization and recombination has the effect of producing charge
states up to Xe8+ in substantial quantities, and of nearly doubling the energy per atom
absorbed by the cluster.

5.4

Cluster Dynamics during the Laser Pulse
As the cluster absorbs energy from the laser field, some of the electrons become

so energetic that they are no longer bound to the cluster. In addition, hydrostatic forcesfrom the hot electrons and Coulomb repulsion in the now positively charged cluster cause
the cluster to expand and cool. These in turn affect the microscopic processes inside thecluster, since all such processes depend on the concentrations of charge species within
the cluster. Collisional ionization and recombination are also sensitively dependent on

84
the temperature of the electron gas relative to the electron binding energy.
The assumption that the electron gas reequilibrates rapidly compared to the rateof photoabsorption can be tested by calculating both rates at all times during the pulse.
The characteristic time scale for energy relaxation in a plasma is controlled by the time
taken for an electron at the plasma temperature to undergo a large-angle deflection asth result of scatterings with other electrons in the plasma. Eliezer, chapter 9.3 [13] gives
this characteristic time as (using atomic units, and putting T in units of energy)√
3 6 T 3/2
τ=.
8 πne ln Λ

(5.31)

This thermalization was compared with the characteristic time for interaction
with the laser beam, defined as the inverse of the total rate of stimulated plus inverse
bremsstrahlung processes per electron in the plasma. The resulting ratio, shown in
Figure 5.4 is highest when the photon flux is low, and falls to a minimum of about 10
when the laser pulse reaches its maximum. Thus, the plasma electrons reequilibratemuch faster than they interact with the laser at all times during the cluster evolution.
The cluster’s expansion during the pulse was tracked using a simple model [10]which tracks only the radius of the cluster, the evaporation of electrons away from thecluster, and the loss of heat from the electron gas resulting from both processes. Grossmovement of electrons or spatial inhomogeneity of charge species within the cluster
were not considered, although a prior study [52] had suggested that such processes
may account for the formation of the highly charged ions detected at the Hamburg
experiment.
Since all processes considered in this model are at most linear with respect toion concentration and are integrated over a cluster of finite volume, such quantities as
the total rate of energy absorption or the total rates of Xe n+ ionizing to form Xe(n+1)+
and recombining to form Xe(n−1)+ are proportional to the number of ions of the given
charge state found in the cluster, but do not change in the event that the charge state

85

Inverse Bremsstrahlung Time / Thermalization Time

1e+06

1e+05

10000

1000

100

10

1

-200

-100
0
100
Time (fs) Relative to Pulse Maximum

200

Figure 5.4: Comparison of thermalization time with time to undergo stimulated orinverse bremsstrahlung for an electron in the cluster plasma, calculated for a 1500 atomcluster exposed to a 7 × 1013 W/cm2 , 100 fs pulse, using the Wigner-Seitz cutoff model
for screening. The electrons thermalize much faster than they interact with the laser at
all times during the laser pulse. Figure taken from [64].

86
distribution is inhomogeneous. Thus, at the level of approximation used, the model
would give identical results if it allowed for spatial inhomogeneity of the various charge
species.
Spatial inhomogeneity of electrons, which was not considered in this study, ispotentially very significant. Volumes with higher electron densities should see higher
rates of inverse bremsstrahlung, collisional ionization and recombination. Subsequentstudies [69, 70] performed using a statistical Boltzmann approach [68] found inhomogeneous distributions of both ions and electrons, with a neutral core filled with highly
charged ions and a high density of electrons. Surrounding this core is a shell of lesshighly charged xenon ions. This shell has a net positive charge due to a lower density ofelectrons. Surrounding the cluster is a cloud of electrons which have enough energy toescape the cluster. The combination of high electron density and highly charged ions at
the center of the cluster would likely increase the total energy absorbed by the clustercompared to the current model.
The equation for the radius of the cluster is given by∂2r
Pe + PCoul 1=3,∂t2
nXe mXe r(5.32)

where Pe = ne Te is the electron pressure and PCoul = Q2 /(8πr 4 ) is the Coulomb pressure
resulting from the charge built up as electrons evaporate away from the cluster.
This model of the laser-cluster dynamics also distinguishes between inner andouter ionization. Inner ionization, which takes place due to photoionization and collisional ionization, is the process by which electrons become liberated from their parent
ion and join the cluster plasma, where they can undergo inverse bremsstrahlung heating or collisional ionization or recombination. Outer ionization is the process by whichelectrons with sufficient energy escape the cluster and cease to interact with it.
The rate of evaporation from a Maxwell distribution of electrons can be calculated
knowing the size of the cluster, the mean free path of electrons in the cluster, and the87temperature of the electron plasma. The rate at which electrons escape from the cluster
is then given by
Wf s =

Z ∞vesc

where vesc =

dv

π λe
(12r 2 − λ2e )vf (v)4 r(5.33)p2(Q + 1)/r is the velocity required for an electron to escape from acluster of charge Q,

v2

f (v) = 4πne (2πT )−3/2 v 2 e− 2T

is the Maxwell distribution, and λe is the mean free path in the cluster plasma, given
by
λe =

T2
4πne (Z + 1) ln Λ

for a plasma with average ion charge Z. The Coulomb logarithm, ln Λ, is set equal to
the logarithm of the screening length in our calculation of the mean free path. λ e isconstrained to be no greater than 2r, the diameter of the sphere.
As electrons evaporate from the cluster, the remaining cluster becomes ever more
highly charged, and a correspondingly lower fraction of the Maxwell distribution has
enough energy to escape from the cluster, thereby choking off the evaporation rate.
It is likely that nearly all high-energy electrons detected in the experiment escapeduring this original period of evaporation. As the cluster expands, the temperature of
the electron plasma falls very quickly as electron thermal energy is converted into ion
kinetic energy, while the energy required to escape the cluster falls only as 1/r.
A recent experiment [30] has for the first time measured the energy spectrum
for electrons emitted from rare gas clusters exposed to intense VUV light. They giveejection spectra for 70 atom xenon clusters exposed to a 4.4 × 10 12 W/cm2 pulse of
VUV light at the same photon energy as the original Hamburg experiment, finding
an electron distribution which decreases approximately exponentially according to I =I0 exp (−Ekin /E0 ), with E0 = 8.9 eV.
The energy spectrum of ejected electrons was calculated by stepping through a88laser pulse using small timesteps. At each timestep, the electron density, mean freepath, cluster radius, and plasma temperature were calculated. These parameters werethen used to calculate the rate at which electrons with energy E = Eesc + Ekin escapedfrom the cluster using Eq. (5.33). Integrated through the timescale of a pulse until theevaporation has stopped, this yields an ejected electron spectrum for a single cluster
exposed to the pulse. Since the clusters are located randomly with respect to the centerof the laser pulse, a further spatial integration over the radial dimension of the pulse,22

assuming a Gaussian laser profile I(r) ∝ e−r /σ from 0 to 3 σ was necessary. The lengthof the interaction region in the Hamburg experiment was comparable to the Rayleighrange for the laser; accordingly, the laser intensity along the direction of propagation
was assumed to be constant. After performing the spatial integration, on average .22
electrons per xenon atom were found to have evaporated from the cluster in this way.
The spectrum of ejection energies for these electrons shown in Fig. 5.5, although notexponential, is nevertheless quite similar to the electron spectrum found in Ref. [30].
The largest discrepancy between this calculated spectrum and the spectrum from[30] occurs at low ejection energy. In addition, the current model of the cluster expansion predicts that the majority of electrons will comprise electron plasmas which
remain bound to the cluster ions and become quite cold during the process of expansion. These electrons – the great bulk of the population – would reach the detector at
low energies and after long delay times, further boosting the spectrum at low energies.
However, Laarmann et al note that for Ekin < 2.5 eV, coinciding with the region oflargest discrepancy, the spectrum cannot be evaluated due to large levels of noise in the
background spectra.
Since electrons faster than about 1 eV are ejected from the cluster during thepulse rather than during the slower process of cluster expansion, the ejected electron
spectrum has the potential to serve as a window into the nature of the laser–cluster
interaction. Accordingly, Fig. 5.6 gives the spectra for 1500 atom clusters exposed to a

89
100 fs, 7 × 1013 W/cm2 pulse, and for 2500 atom clusters exposed to a 50 fs, 2.5 × 1013
W/cm2 pulse.
The calculated spectrum was different for different models of plasma screening.
After spatial averaging, 1500 atom clusters exposed to a 100 fs, 7 × 10 13 W/cm2 pulse
were found to eject 0.22 electrons per atom during this early evaporation period using
the Wigner-Seitz cutoff model for the screening length (see section 5.5 for a discussionof plasma screening). Using the Attard model, 0.07 electrons per atom are evaporated
during this period. For 2500 atom clusters exposed to a 50 fs, 2.5 × 10 13 W/cm2 pulse,
the corresponding numbers are 0.13 electrons per atom for the Wigner-Seitz cutoff modeland 0.02 electrons per atom for the Attard model. In contrast to this, the Hamburg
experiment measured an average charge per ion of 2.98. Hence, the electrons which
comprise these ejected electron spectra correspond to only a few percent of all freeelectrons at the time when the expanding clusters reach the detector.

5.5

Nonideal Plasma Screening
As shown in Fig. 5.8, when plasma screening of the Xe ions becomes strong

enough to allow photoionization of Xe+ into Xe2+ , large numbers of extremely lowenergy electrons are added to the plasma. As a result, the ratio of electron kinetic energy to electrostatic potential energy falls dramatically, the Debye length of the plasma
falls abruptly below the Wigner-Seitz radius of xenon, and the plasma enters a regimeof strong correlation. In this regime, a number of the assumptions of Debye-Hückel
screening model break down, and the Debye length loses its meaning as a screening
distance [18]. If the plasma cools sufficiently, screening lengths can become complex,
and result in oscillatory electron–ion correlation functions [31, 3].
Another possibly important effect of the strongly coupled plasma was identified ina recent study [25], which has identified electron dynamics in a strongly coupled plasma
as having a very large impact upon rates of many-body recombination and hence upon

90

1

Ejected Electrons (arb. units)

Laarmann et al
Wigner-Seitz Screening
Attard Screening
0.75

0.5

0.25

0

0

10

20
Ejected Electron Energy (eV)30

40

Figure 5.5: Ejected electron spectrum. Comparison between data from [30] andspatially-averaged spectra calculated using 70 atom clusters exposed to a 4.4 × 10 12
W/cm2 , 100 fs pulse for two different models of plasma screening. The Wigner Seitz
cutoff model uses the ordinary Debye length as the screening radius, but the screening radius is not allowed to fall below xenon’s Wigner-Seitz radius at liquid density,
4.64 bohr. The Attard model of screening calculates the screening radius according
to equation (5.34), discussed in Section 5.5. The spectrum calculated using xenon’s
Wigner-Seitz radius as a minimum screening distance displays a strong similarity to theexperimental curve. The intensity of the experimental spectra is arbitrary; magnitudes
were chosen by setting each curve equal at the beginning of the exponential tail in theexperiment. Figure taken from [64].

91

1
a)

Wigner-Seitz Screening
Attard Screening

Ejected Electrons (arbitrary units)

0.8

0.6

0.4

0.2

0

0

20

40
60
Kinetic Energy (eV)

80

100

1
b)

Wigner-Seitz Screening
Attard Screening

Ejected Electrons (arbitrary units)

0.8

0.6

0.4

0.2

0

0

20

40
60
Kinetic Energy (eV)

80

100

Figure 5.6: Ejected electron spectra, calculated for the two sets of parameters andthe two models of screening. The Wigner Seitz cutoff model uses the ordinary Debye
length as the screening radius, but the screening radius is not allowed to fall below
xenon’s Wigner-Seitz radius at liquid density, 4.64 bohr. The Attard model of screening
calculates the screening radius according to equation (5.34), discussed in section 5.5.
a) Nature parameters: 1500 atom clusters exposed to a 100 fs, 7 × 10 13 W/cm2 pulse.
b) Thesis parameters: 2500 atom clusters exposed to a 50 fs, 2.5 × 10 13 W/cm2 pulse.
Since electrons faster than about 1 eV are ejected from the cluster during the pulse
the ejection spectra could serve as a window into the dynamics of the laser-cluster
interaction. Figure taken from [64].

92

30
a)

Total Energy Absorbed
Electron Kinetic Energy

Energy (hartree) per Atom

25

20

15

10

5

0

-200

1

-100
0
100
Time (fs) Relative to Pulse Maximum

0+

b)

1+

0+
1+
2+
3+
4+
5+
6+
7+
8+

Fraction of Population

0.8
2+ 3+
0.6

4+

200

5+

0.4
6+
8+

0.2

0

7+

-200

-100
0
100
Time (fs) Relative to Pulse Maximum

200

Figure 5.7: The effects of collisional ionization and recombination are to allow theformation of charge states beyond Xe3+ Pictured is the time evolution of a single 1500
atom cluster exposed to a 100 fs, 7 × 1013 W/cm2 pulse. These states enhance therate of inverse bremsstrahlung heating. As the plasma expands and cools, the chemical
equilibrium shifts toward lower charge states on a timescale much longer than the laser
pulse, until decreasing plasma density causes recombination and ionization rates to go
to zero. a) Energy absorbed vs. time b) Ionic population vs. time during laser pulse.
Figure taken from [64].

93energy absorption by the cluster as the recombined ions undergo multiple episodes ofphotoionization.
Most calculations performed in this chapter were performed using xenon’s WignerSeitz radius at liquid density as a minimum value below which the screening was notallowed to fall. Clearly, with the precise nature of screening unknown in the strongly correlated regime, the method of calculating atomic properties based on a Debye-screened
atomic potential acquires a corresponding uncertainty. Accordingly, as an attempt toestimate this uncertainty, the evolution of the cluster was calculated using different
models for the screening length.
The simplest approximation applied xenon’s Wigner-Seitz radius at liquid density as a minimum value below which the screening length was not allowed to fall. A
second model, proposed by Attard [3], deals with ions having a nonzero radius. Strictly
speaking, the Debye-Hückel model for plasma screening is invalid except in the limit
of ions which have zero size. Attard has shown that in the case where ions have anonzero hard-sphere radius d, the screening length λ = 1/κ differs from the classical
Debye-Hückel length λD = 1/κD according to
κD
κ= p
.
2
1 − (κD d) /2 + (κD d)3 /6(5.34)

This effect becomes important in the domain where λD ≤ d.
Qualitatively, the effect of considering screening lengths in this model which areshorter than the Wigner-Seitz radius is twofold. First, the tighter screening slightly
decreases inverse bremsstrahlung heating. Secondly, it allows photoionization of Xe 3+
and higher charge states. Directly substituting the Attard screening length for the
Debye length with Wigner-Seitz cutoff therefore gives some insight as to how sensitive
the results in this chapter are to different models of the ionic potential under verystrong screening. As can be seen in Fig. 5.9 the Attard screening model has a relatively
small impact on the prediction for the energy absorbed by the cluster. More prominent

94

20
Wigner-Seitz Cutoff
Attard Screening

Screening Length (bohr)

15

10

5

0

-200

-100
0
100
Time (fs) Relative to Pulse Maximum

200

Figure 5.8: The interaction of plasma screening with atomic potentials is unknown as
the screening length becomes very short. Here the screening length vs time is given for
two simulations of a 1500 atom cluster exposed to a 100 fs, 7 × 10 13 W/cm2 pulse, using
two models for screening. In the first model, the screening length is not allowed to fall
below xenon’s Wigner-Seitz radius at liquid density. The second model for screening
uses a formula given by Attard. Figure taken from [64].

95

30
a)

Total Energy Absorbed
Electron Kinetic Energy

Energy (hartree) per Atom

25

20

15

10

5

0

-200

1

-100
0
100
Time (fs) Relative to Pulse Maximum

0+

b)

1+

0.8
Fraction of Population

200

8+

0+
1+
2+
3+
4+
5+
6+
7+
8+

0.6

2+

0.4

3+ 4+ 5+
6+

0.2
7+
0

-200

-100
0
100
Time (fs) Relative to Pulse Maximum

200

Figure 5.9: Near the center of the pulse, calculating the evolution of the cluster using
Attard screening yields screening lengths shorter than the Wigner-Seitz cutoff, allowing
easier formation of high charge states than when the evolution is calculated for WignerSeitz screening, shown in Figure 5.7. For a 1500 atom cluster exposed to a 7 × 10 13
W/cm2 , 100 fs pulse: a) Energy absorbed vs time for the Attard screening model, b)
Charge species population vs time for the Attard screening model. Figure taken from
[64].

96

80
a)

Electron-Electron Coupling
Electron-Ion Coupling

Coupling Parameter

60

40

20

0

-200

Coupling Parameter

4000

b)

-100
0
100
Time (fs) Relative to Pulse Maximum

200

Electron-Electron Coupling
Electron-Ion Coupling

3000

2000

1000

0

-200

-100
0
100
Time (fs) Relative to Pulse Maximum

200

Figure 5.10: Plasma coupling parameter vs. time for the two models of plasma screening.
3/2
1
and Γei = ZΓee where the average
The coupling parameters are defined by Γee = aT3
)1/3 and Z is the average charge
distance between electrons a is given by a = ( 4πneof the ions. The plasma becomes very strongly coupled early in the pulse, but thestrength of the coupling decreases as the plasma absorbs more energy in the course of
the cluster heating. a) Coupling parameters vs. time using Wigner-Seitz cutoff. b)
Coupling parameters vs. time for the Attard screening model. Figure taken from [64].

97
is the formation of higher charge states, which is abetted by the reduced ionization
potentials resulting from the tighter screening in the Attard model. Figure 5.10 shows
the plasma coupling parameter, a measure for the nonideality of a plasma, for the twomodels, demonstrating that the Attard screening model gives rise to a more strongly
coupled plasma than the pure Debye model. In addition, the two models give different
populations for the various charge states at the end of the pulse; however, the combined
effects of the cluster expansion and spatial averaging over the beam profile act to destroy
much of this information.
As a further test of our model’s sensitivity to plasma screening lengths, the evolution of the cluster was calculated in the limit of weak screening. To do this, thescreening length in the cluster was fixed at 100 bohr for the entire duration of the laser
pulse. In this limit, ionization potentials are unchanged from their values in the absenceof screening, and photoionization beyond Xe+ is impossible. In this model, formation ofstates with charge 2 or higher must come entirely from collisional ionization. The resultsof this constraint can be seen in Figure 5.13. Neglecting the effects of plasma screening
in this way inhibits the formation of high charge states in the cluster, yielding virtuallyno Xe7+ or Xe8+ . After the initial photoionization of neutral xenon, the plasma undergoes a period of slow heating while the laser intensity builds. Near the maximum of
the pulse, the plasma becomes energetic enough to ionize the higher charge states with
their unscreened ionization potentials, and the rate of inverse bremsstrahlung heating
increases rapidly. The total energy absorbed falls from 25 hartree per atom in the caseof the Wigner-Seitz cutoff model to 11.5 hartree in the limit of no screening.

5.6

Hydrogenic Model of Inverse Bremsstrahlung
Most previous approaches to the problem of laser-cluster interactions have consid-

ered the ionic potential seen by the electron as a pure Coulomb potential. This is not an
unreasonable approximation: as the charge of the ion increases, the difference between

98
inverse bremsstrahlung cross sections calculated using Herman Skillman potentials andcross sections calculated using Coulomb potentials is much smaller than in the case of
the bare ion. This can be seen in figure 5.11, which contrasts inverse bremsstrahlungcross sections calculated using Coulomb and Herman-Skillman potentials for ions ofcharge 5.
As can be seen in figure 5.7b, when the laser reaches maximum intensity, mostof the cluster has been ionized to such high charge states. Thus, models of the inverse
bremsstrahlung process which use Coulombic potentials should be able to see comparable levels of heating to those using cross sections derived using Herman-Skillman
potentials.
To investigate this proposition, the laser-cluster interaction for a 1500 atom cluster exposed to a 100 fs, 7 × 1013 W/cm2 was simulated using the model presented inthis chapter, but with a physical picture chosen to emulate that of Siedschlag and Rost[52]. The simulation used inverse bremsstrahlung cross sections calculated with Debyescreened Coulomb potentials. Ionization potentials and photoionization cross sections
were unchanged from the other simulations. The unaltered Debye length was used as
the screening length. Collisional ionization and recombination were not considered.
The results of this simulation are presented in figure 5.12. Levels of energy absorption were found which were very comparable to those in the model using HermanSkillman potentials, but the behavior of ionic populations with time was very different.
Xe7+ and Xe8+ , which make up almost half of the population of the cluster at the endof the pulse in the Herman-Skillman model, were present in negligible quantities.
Both differences between the two physical pictures are attributable to the effectsof collisional ionization and recombination. Recombination slows the growth of highcharge state populations by allowing some photoionized ions to recombine into a lower
charge state, while collisional ionization allows the population of charge states whichcannot be created via sequential photoionization.

Inverse Bremsstrahlung Cross Section (Mb/hartree)

99

Coulomb Potential
Herman-Skillman Potential
10000

1000

0

5
Initial Electron Energy (hartree)

10

Figure 5.11: Inverse bremsstrahlung cross sections [Eqs. 5.13 and 5.14] calculated for
an electron in the field of a purely Coulombic 5+ potential and for an electron in thefield of a Xe Herman-Skillman atomic potential of the same charge. In comparison with
figure 5.1, it can be seen that at higher charge states, the impact of atomic structure
on inverse bremsstrahlung cross sections is decreased. Figure taken from [64].

100
20

a)
Total Energy Absorbed
Electron Kinetic Energy

Energy (hartree) per Atom

15

10

5

0

-200

1

-100
0
100
Time (fs) Relative to Pulse Maximum

0+

b)

6+
0+
1+
2+
3+
4+
5+
6+
7+
8+

0.8
Fraction of Population

200

0.6

0.4

1+

0.2

2+
3+
0

5+
4+

5+

-200

-100
0
100
Time (fs) Relative to Pulse Maximum

200

-200

-100
0
100
Time (fs) Relative to Pulse Maximum

200

20
c)

Screening Length (bohr)

15

10

5

0

Figure 5.12: Simulation of the laser-cluster interaction using a physical model takenfrom [52]. In this model, inverse bremsstrahlung cross sections are calculated usinghydrogenic potentials and all high charge states are produced via sequential photoionization. Collisional ionization and recombination are not considered. a)Energy absorbed
vs time. b) Charge state population vs time. c) Debye length vs time. Figure takenfrom [64].

101

15

Energy (hartree) per Atoma)Total Energy Absorbed
Electron Kinetic Energy

10

5

0

-200

-100
0
100
Time (fs) relative to Pulse Maximum

200

1
0+

b)

Fraction of Population

0.8

1+

0+
1+
2+
3+
4+
5+
6+

2+

0.6
3+
4+

0.4

5+

0.2

6+
0

-200

-100
0
100
Time (fs) Relative to Pulse Maximum

200

Figure 5.13: Evolution of the cluster calculated in the limit of long screening length,
λD = 20 bohr for a 1500 atom cluster exposed to a 100 fs, 7 × 10 13 W/cm2 pulse. Weakscreening prevents ionization potential lowering due to plasma effects, and precludes
photoionization past Xe+ . The effect of neglecting plasma screening effects is to reduceformation of high charge states and to reduce the total energy absorbed by the cluster.
The pulse and cluster parameters are identical to those used in Figures 5.3, 5.7 and 5.9.
a) Energy absorbed vs. time b) Ionic population vs. time during laser pulse. Figure
taken from [64].

102

5.7

Summary and Conclusions
When a xenon cluster is irradiated by intense VUV light, there are four phases in

its evolution. In the first phase, electrons are liberated from the xenon atoms and forma plasma. As the number of free electrons grows, the screening length of the plasma
shrinks.
Once the screening length of the plasma reaches 10.6 bohr, Xe1+ can undergophotoionization into Xe2+ . This results in the addition of large numbers of low-energy
electrons to the plasma, cooling it and decreasing the screening length still further. The
ratio of kinetic energy to potential energy falls dramatically, and the plasma temporarily becomes strongly coupled. Ionization potentials for higher charge states fall withincreased screening, facilitating their creation.
In the third phase, the plasma undergoes rapid inverse bremsstrahlung heating.
High charge states are formed through collisional ionization and recombination, and thecluster becomes charged as energetic electrons evaporate away from its surface. The
charge state distribution shifts rapidly toward higher charges, with the average ioniccharge reaching 5.5 at the pulse peak. This distribution changes only slowly on thetimescale of the pulse.
Finally, the cluster expands due to the pressure of the electron gas and the cluster’s own charge. As the cluster expands, the electron plasma cools and becomes morediffuse. Screening lengths increase, and charge state equilibrium shifts toward lower
charge states.
Of these four phases, the model presented here describes the first and third phases
well; the second more crudely. The dynamics of the expanding cluster are a challenging
problem in their own right, and demand a treatment more sophisticated than the current
simple homogeneous expansion model.
For strongly coupled plasmas, it is unclear whether this treatment of plasma

103screening adequately describes the potential seen by scattering or photoionizing electrons. As the Debye length falls below the Wigner-Seitz radius, the interaction of
screening effects due to inner-shell electrons and effects due to screening by continuum
plasma electrons should be considered. It is known that the screening length diverges
from the Debye length in this limit, but the precise nature of the electron–ion potential
is unknown.
There is some difficulty in comparing these results to the Hamburg experiment,
due to experimental uncertainty in laser intensity, temporal profile, spatial profile, andcluster size. Whereas in the Nature paper the Hamburg group described the laser pulse
as 100 fs, 7 × 1013 W/cm2 incident on 1500 atom xenon clusters, Wabnitz’s thesis [62]subsequently describes these pulses as 50 fs, 2.5 × 1013 W/cm2 pulses incident on 2500
atom clusters. In addition, the temporal profile of the laser pulses is not Gaussian,and varies in an unpredictable way from pulse to pulse due to the nature of the SASE
amplification process, which starts from shot noise.
This model also has difficulty explaining the properties of the clusters long after
the laser-cluster interaction is over. As the clusters expand and cool, they continue
to undergo collisional ionization and recombination. The distribution of charge states
measured at the experimental detectors bears no simple relationship to the distribution
calculated at the end of the pulse. The homogeneous model of the cluster expansion
implicitly requires that all charge states in the same cluster have the same average kinetic
energy; this obviously conflicts with the quadratic dependence of energy vs charge state
detected in the Hamburg experiment. Also, it is likely that high charge states escapethe cluster more quickly than low charge states, spending less time in regions of highelectron density and having less opportunity to recombine. Thus, a more sophisticatedmodel of the cluster expansion is necessary in order to predict final charge state andionic energy distributions with confidence for comparison with experiment.
At the center of a gaussian laser pulse using parameters taken from the Nature

104
paper and a Wigner-Seitz debye length cutoff, each cluster absorbs on average 682 eV
per atom. At a distance of 3 sigma from the center of such a gaussian pulse, eachcluster absorbs only .4 eV per atom. Spatial averaging over the gaussian pulse profile
from 0 to 3 sigma gives an average of 195 eV per atom absorbed. Using the less intenseparameters from Wabnitz’s thesis gives 219 eV per atom at the center, 0.2 eV per atom
at 3 sigma, and 65 eV per atom on spatial averaging.
Using a time of flight detector which could detect only charged ions, Wabnitz et
al. reported an average ion energy of 400 eV, subsequently revised to 650 eV.
Clearly, a spatial average such as the one performed here could be altered byaveraging over a different beam profile or by changing the limits of the radial averageand including more clusters which are exposed to only a tiny fraction of the beam’s
peak intensity. It is also clear that most of the atoms in the clusters which are exposedto very small fractions of the peak intensity will never be ionized and thus would notregister in a time-of-flight ion detector such as was used in the Hamburg experiment.
Thus, in the absence of better information about the beam’s spatial and temporal profile
and a more comprehensive model of the cluster expansion after the conclusion of thelaser pulse, it is impossible to make precise comparisons between this model and theHamburg results.
Nevertheless, this model of the laser-cluster interaction explains some surprisingfeatures of the laser-cluster interaction in the VUV regime quite well. Primary among
these is the surprising efficiency by which the clusters absorb photons. Second, the highcharge states observed in the Hamburg experiment here emerge quite naturally. Third,
the same model is able to calculate the early electron ejection spectrum measured in [30]and achieved great similarity to experiment, despite a cluster size and pulse intensity
which differ significantly from those of the original Hamburg experiment. Such spectra
can depend strongly on the model of plasma screening or the precise parameters of
the experiment, and can therefore serve as a possible window into the nature of the105
laser-cluster dynamics during the time period of the pulse.

Chapter 6

Summary, Conclusions and Outlook

This thesis has investigated the ways in which atomic and molecular physics
are affected by the presence of an intense laser field. Chapter 2 examined how thefamiliar concept of a shape resonance is altered by the intense and time-varying fieldof the driving laser, while chapter 3 investigated how molecular scattering states affectexperimental observables such as the high harmonic spectrum, as well as the implications
for attempts to reconstruct molecular properties from such information. Chapter 4
investigated how the internal degrees of freedom in a molecule can measurably affecthigh harmonic experiments. Chapter 5 examined the physically different scenario of alaser-cluster interaction, in which atomic properties are altered by interaction with acluster plasma rather the driving laser.
As is not surprising for such a young field, the work presented in this thesis
can hardly be said to encompass everything which could be done. The most obvious
omission from this thesis is a satisfactory treatment of molecular ionization in a strong
field. Effects such as polarization of the molecular core or ionization from multiple
orbitals have barely been considered here, and in nowhere near the depth they deserve.
In the longer run, more elaborate treatments of the molecular wavefunction offer
great potential as ways to look at molecular dynamics on an ultrafast timescale. A simple variation of the quantum interference experiment described in chapter 4 would be
to use the first laser pulse to start some chemically interesting process such as a change

107
of molecular conformation and then to use the second pulse as an interferometric probeof the vibrational wavefunction dynamics. Another particularly attractive possibility
would be to investigate vibronic dynamics in the vicinity of a conical intersection byionizing a symmetric molecule near the symmetry point and examining how the vibrational wavefunction affected the electronic dynamics in the ion. Such experiments could
close the loop: atomic and molecular ideas, once understood in the exotic regime ofintense laser physics, could in turn provide a new window to look at chemical dynamics
on previously unreachable timescales.
To return to the train metaphor from the introduction, a bright glare shining inthrough the windows does not fundamentally alter the interactions between the passengers, but it may serve to illuminate those interactions better than a static seating
arrangement could. The fact that you avoid the boor to your front and gravitate toward
your friends in the back could be an accident in any particular seating chart, but would
be rapidly apparent after the cabin has rearranged itself a few times. Likewise, the goalof dealing with atoms and molecules in intense laser fields is not merely to do atomic
and molecular physics in an exotic regime, but also as a means of learning about normal
molecules in normal environments.

Bibliography[1] N. Wagner, private communication.
[2] In practice, matrix elements involving large values of l were found to converge very
quickly to zero. Because of this, the total transition rate from states of energy E
to states of energy E 0 , and hence the heating, is unaffected by the existence of acutoff.
[3] Phil Attard. Phys. Rev. E, 48:3604, 1993.
[4] M. Aymar, C. H. Greene, and E. Luc-Koenig. Rev. Mod. Phys., 68:1015, 1996.
[5] I.B. Bersuker. The Jahn-Teller effect and vibronic interactions in modern chemistry.
Plenum Press, 1984.
[6] V. Carravetta, Y. Luo, and H. Agren. Chem. Phys., 174:141, 1993.
[7] P. B. Corkum. Phys. Rev. Lett., 71:1994, 1993.
[8] C. Deiss, N. Rohringer, J. Burgdoerfer, E. Lamour, C. Prigent, and D. Vernhet.
Journal de Physique IV, 138:55, 2006.
[9] T. Ditmire, T. Donnelly, R. W. Falcone, and M. D. Perry. Phys. Rev. Lett., 75:3122,
1995.
[10] T. Ditmire, T. Donnelly, A. M. Rubenchik, R. W. Falcone, and M. D. Perry. Phys.
Rev. A, 53:3379, 1996.
[11] T.H. Dunning Jr. Gaussian basis sets for use in correlated molecular calculations.
i. the atoms boron through neon and hydrogen. Journal of Chemical Physics,
90(2):1007–1023, 1989.
[12] A. R. Edmonds. Angular Momentum in Quantum Mechanics. Princeton University
Press, Princeton, New Jersey, 1996.
[13] Shalom Eliezer. The Interaction of High-Power Lasers with Plasmas. Institute ofPhysics Publishing, Bristol and Philadelphia, 2002.
[14] B. D. Esry and Chris H. Greene. unpublished, 1995.
[15] S. Estreicher and T. L. Estle. Phys. Rev. B, 31:5616, 1985.

109[16] M. J. Frisch et al. Gaussian 98. Gaussian Inc., Pittsburgh, PA, 1998.
[17] V. Ayvazyan et al. Phys. Rev. Lett., 88:104802, 2002.
[18] Vladimir E. Fortov and Igor T. Iakubov. The Physics of Non-Ideal Plasma. World
Scientific, Singapore, 2000.
[19] Harald Friedrich. Theoretical Atomic Physics. Springer Verlag, Berlin, 1991.
[20] H. Goldstein, C. Poole, and J. Safko. Classical Mechanics, Third Edition. Addison
Wesley, 2002.
[21] M. C. Gutzwiller. Chaos in Classical and Quantum Mechanics. Springer, New
York, 1990.
[22] F. Herman and S. Skillman. Atomic Structure Calculations. Prentice-Hall, Englewood Cliffs, N.J., 1963.
[23] M. H. R. Hutchinson, T. Ditmire, E. Springate, J. W. G. Tisch, Y. L. Shao, M. B.
Mason, N. Hay, and J. P. Marangos. Phil. Trans. R. Soc. Lond. A, 356:297, 1998.
[24] J. Itatani, J. Levesque, D. Zeidler, Hiromichi Niikura, H. Pepin, J. C. Kieffer, P. B.
Corkum, and D. M. Villeneuve. Nature, 432:867, 2004.
[25] C. Jungreuthmayer, L. Ramunno, J. Zanghellini, and T. Brabec. J. Phys. B.,
38:3029, 2005.
[26] S. Kazamias and P. Balcou. Physical Review A, 69:063416, 2004.
[27] D. A. L. Kilcoyne, S. Nordholm, and N. S. Hush. Chemical Physics, 107:197, 1986.
[28] J. L. Krause, K. J. Schafer, and K. C. Kulander. Phys. Rev. A, 45:4998, 1992.
[29] W. L. Kruer. The Physics of Laser Plasma Interactions. Westview Press, Boulder,Colorado, 2003.
[30] T. Laarmann, M. Rusek, H. Wabnitz, J. Schulz, A. R. B. de Castro, P. Guertler,
W. Laasch, and T. Moeller. Phys. Rev. Lett., 95:063402, 2005.
[31] Benjamin P. Lee and Michael E. Fisher. Phys. Rev. Lett., 76:2906, 1996.
[32] M. Lein. J. Phys. B, 40:R135, 2007.
[33] M. Lewenstein, P. Balcou, M.Y. Ivanov, A. L’Huillier, and PB Corkum. Phys. Rev.
A, 49:2117, 1994.
[34] M. Lewenstein, Ph. Balcou, M. Yu. Ivanov, A.L’Huillier, and P. B. Corkum. Phys.
Rev. A, 49:2117, 1994.
[35] M. Lezius, S. Dobosz, D. Normand, and M. Schmidt. Phys. Rev. Lett., 80:261,
1998.
[36] Wolfgang Lotz. Zeitschrift fur Physik, 216:241–247, 1968.

110[37] W. H. Miller. J. Phys. Chem. A, 105:2942, 2001.
[38] W. Moffitt and W. Thorson. Phys. Rev., 108:1251, 1957.
[39] T. Morishita, A. T. Le, Z. Chen, and C. D. Lin. preprint.
[40] H. Nakamura. J. Theor. Comp. Chem., 4:127, 2005.
[41] L. Nugent-Glandorf, M. Scheer, D. A. Samuels, V. M. Bierbaum, and S. R. Leone.
J. Chem. Phys., 117:6108, 2002.
[42] S. Patchkovskii, Z. Zhao, T. Brabec, and D. M. Villeneuve. J. Chem. Phys.,
126:114306, 2007.
[43] J. Plenge, C. Nicolas, A. G. Caster, M. Ahmed, and S. R. Leone. J. Chem. Phys.,
125:133315, 2006.
[44] William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical Recipes in Fortran 77. Cambridge University Press, Cambridge,
England, 2001.
[45] T. N. Rescigno and C. W. McCurdy. Physical Review A, 62:032706, 2000.
[46] Manuel Rotenberg. The 3-j and 6-j Symbols. Technology Press, Cambridge, Massachusetts, 1959.
[47] R. Santra and A. Gordon. Phys. Rev. Lett., 96:073906, 2006.
[48] R. Santra and C. H. Greene. Phys. Rev. A, 70:053401, 2004.
[49] Robin Santra and Chris H. Greene. Phys. Rev. Lett., 91:233401, 2003.
[50] A. Scrinzi and N. Elander. J. Chem. Phys, 98:3866, 1992.
[51] Y. L. Shao, T. Ditmire, J. W. G. Tisch, E. Springate, J. P. Marangos, and M. H. R.
Hutchinson. Phys. Rev. Lett., 77:3343, 1996.
[52] Christian Siedschlag and Jan-Michael Rost. Phys. Rev. Lett., 93:043402, 2004.
[53] O. Smirnova, M. Spanner, and M. Ivanov. J. Phys. B, 39:S307, 2006.
[54] G. Snell, U. Hergenhahn, N. Mueller, M. Drescher, J. Viefhaus, U. Becker, and
U. Heinzmann. Phys. Rev. A, 63:032712, 2001.
[55] A. F. Starace. Phys. Rev. B, 5:1773, 1972.
[56] J. Sugar. Physical Review B, 5(5):1785, 1972.
[57] S. Tonzani. Comp. Phys. Comm., 176:146, 2007.
[58] S. Tonzani and C. H. Greene. J. Chem. Phys., 122:014111, 2005.
[59] S. Tonzani and C. H. Greene. J. Chem. Phys., 124:054312, 2006.

111[60] H. Wabnitz, L. Bittner, A. R. B. de Castro, R. Doehrmann, P. Guertler, T. Laarmann, W. Laasch, J. Schulz, A. Swiderski, K. von Haeften, T. Moeller, B. Faatz,
A. Fateev, J. Feldhaus, C. Gerth, U. Hahn, E. Saldin, E. Schneidmiller, K. Sytchev,
K. Tiedtke, R. Treusch, and M. Yurkov. Nature, 420:482, 2002.
[61] H. Wabnitz, A. R. B. de Castro, P. Guertler, T. Laarmann, W. Laasch, J. Schulz,
and T. Moeller. Phys. Rev. Lett., 94:023001, 2005.
[62] Hubertus Wabnitz. Interaction of Intense VUV Radiation from a Free-Electron
Laser with Rare Gas Atoms and Clusters. PhD thesis, Universitat Hamburg, 2003.
[63] N. L. Wagner, A. Wuest, I. P. Christov, T. Popmintchev, X. Zhou, M. M. Murnane,and H. C. Kapteyn. Proc. Natl. Acad. Sci. U.S.A., 103(36):13279, 2006.
[64] Z. Walters, R. Santra, and C. H. Greene. Phys. Rev. A, 74:043204, 2006.
[65] Z. Walters, S. Tonzani, and C. H. Greene. J. Phys. B, 40:F277, 2007.
[66] Z. B. Walters, S. Tonzani, and C. H. Greene. J. Phys. Chem. A, 112:9439–9447,
2008.
[67] Ya. B. Zel’dovich and Yu. P. Raizer. Physics of Shock Waves and High-Temperature
Hydrodynamic Phenomena. Academic Press, New York and London, 1966.
[68] B. Ziaja, A. R. B. de Castro, E. Weckert, and T. Moeller. European Physics Journal
D, 40:465, 2006.
[69] B. Ziaja, H. Wabnitz, E. Weckert, and T. Moeller.
10:043003, 2008.

New Journal of Physics,[70] B. Ziaja, H. Wabnitz, E. Weckert, and T. Moeller. European Physics Letters,
82:24002, 2008.