ULTRAFAST SCIENCE: A MULTISCALE MODELING APPROACH TO
FEMTOSECOND ELECTRON DIFFRACTION AND ITS APPLICATIONS
By
Jenni Minttu Eleonora Portman

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Physics - Doctor of Philosophy
2014

ABSTRACT
ULTRAFAST SCIENCE: A MULTISCALE MODELING APPROACH TO
FEMTOSECOND ELECTRON DIFFRACTION AND ITS APPLICATIONS
By
Jenni Minttu Eleonora Portman
The focus of this work is the study of processes at the edge of the current space and time
resolutions. This includes both efforts in the development of an ultrafast electron microscope
(UEM) and in the study of correlated electron systems that reflect measurements taken with
this instrument.
The development of a reliable ultrafast electron diffraction and imaging system requires a
low emittance source of photoemitted electrons and an understanding of how the properties
of the generated bunch depend on the photocathode properties. In order to gain more
understanding of this process, we combine the so-called three-step photoemission model with
N-particle electron simulations. By using the Fast Multipole Method to treat space charge
effects, we are able to follow the time evolution of pulses containing over 106 electrons and
investigate the role of laser fluence and extraction field on the total number of electrons that
escape the surface as well as virtual cathode physics and the limits to spatio-temporal and
spectroscopic resolution originating from the image charge on the surface and from the profile
of the exciting laser pulse. The results of these simulations are compared to experimental
images of the photoemission process collected using the shadow imaging technique. By
contrasting the effect of varying surface properties (leading to expanding or pinned image
charge) and laser profiles (Gaussian, uniform and elliptical) under different extraction field
strengths and numbers of generated electrons, we quantify the effect of these experimental
parameters on macroscopic pulse properties such as emittance, brightness (4D and 6D),

coherence length and energy spread. Based on our results, we outline optimal conditions of
pulse generation for UEM systems.
With our knowledge of the photoemitted pulse properties, we also present our development of a design for the whole UEM column using the Analytic Gaussian model. We
summarize the derivation of the equations governing this mean field model and show how
the contributions due to the photoemission gun and the relativistic motion of the electrons
can be added to this formalism to make it applicable for our system. We then explain the
procedure used for optimizing the lens and RF cavity strengths and analyze both the effect of each separate optical element and their role in the column. We discuss the limits
of this model and calculate the achievable temporal and spatial resolutions under different
photoemission conditions.
We conclude the present work with an investigation of Tantalum Disulphide (TaS2 ), a
material that presents interesting ultrafast phenomena that can be probed using the UEM.
TaS2 is a transition-metal layered compound that for T <190 K displays a commensurate
charge density wave (C-CDW) phase characterized by insulating behavior with the opening of
a gap at the Fermi energy. To better understand the C-CDW phase and explain its quantum
mechanical origin, we perform density functional theory calculations of the electronic band
structure of 1T-TaS2 and quantify the effect that spin orbit coupling and Hubbard repulsion
have on the ground state. Our results show that neither of these interactions is sufficient to
reproduce the insulating gap seen in experiment, an observation which is confirmed by our
calculation of the phonon band structure and absorption spectrum. We also consider the
effect of different stacking configurations of the TaS2 layers and find evidence of gap opening
for bilayers in the presence of disordered stacking.

Yet all experience is an arch wherethro’
Gleams that untravell’d world whose margin fades
For ever and forever when I move.
— Tennyson, Ulysses

iv

ACKNOWLEDGMENTS

Science is the fruit of a collaborative effort, so I would like to recognize the people who have
shaped me as a researcher and impacted my work in meaningful ways.
First, I would like to thank my advisor Phil for always having time and for giving me
freedom in my work to pursue and grow my own scientific interests. I would like to thank
Chong Yu and Martin, it has been very valuable for me to be part of our collaboration and
be able to learn from you. I also thank the other members of my committee, Remco and
Alex for taking the time to understand my work. To the other students working on the
UEM project, He Zhang, Zhensheng Tao, and Kiseok Chang, it has been a pleasure to work
with you. I would also like to recognize the help and advice that both Prof. Mahanti and
Dat Thanh Do have given me. As all my simulations were carried on the HPCC clusters at
MSU, I cannot help but mention them and thank their staff, in particular Dirk Colbry.
Next is a big thank you to my friends and the other members of the group. Dan, Connor,
Jim and Shannon your friendship is very important for me. Thank you also to Gift, Bin and
Xukun for making our office a wonderful place to work in.
I also want to mention my family, my mother Anneli, who beat me in getting a PhD by a
few months, my father Kim, who let me mess around with computers since an early age, my
siblings Joonatan, Sonja, Kristofer and my sister-in-law Valentina for being awesome. Last
but not least, my husband Giuseppe, who is my bearded number one fan. He has always
encouraged my work, given me excellent insight and advice, and supported me in many
different ways. I couldn’t have asked for a better partner for my journey.
I acknowledge support from US National Science Foundation under Grant No. NSF-DMR
1126343 for the work included in this dissertation.
v

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 1 Introduction . . . . . . . . . . .
1.1 Why study ultrafast science . . . . . .
1.2 Pump probe microscopy . . . . . . . .
1.3 Ultrafast Electron Microscope at MSU
1.4 Overview . . . . . . . . . . . . . . . . .

.
.
.
.
.

1
1
2
4
5

Chapter 2 Generation of the electron pulse . . . . . . . . . . . . . . . .
2.1 Importance of the initial stages of pulse generation . . . . . . . . . . .
2.2 Modeling the photoemission process . . . . . . . . . . . . . . . . . . . .
2.3 Details of the simulations . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Treatment of the positive surface field . . . . . . . . . . . . . . . . . . .
2.4.1 Equations for the surface field: Gaussian, Uniform and Elliptical
2.4.2 Time evolution of the positive field . . . . . . . . . . . . . . . .
2.5 Pulse time evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6 Optimal pulse generation . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.1 Gaussian pulse . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.2 Effect of image charge pinning . . . . . . . . . . . . . . . . . . .
2.6.3 Effect of transverse laser profile . . . . . . . . . . . . . . . . . .
2.7 Guidelines for generation of high brightness electron beams . . . . . . .

. . .
. . .
. . .
. . .
. . .
cases
. . .
. . .
. . .
. . .
. . .
. . .
. . .

7
7
9
14
17
17
20
22
28
29
31
33
35

Chapter 3 Analytic Gaussian Model . . . . . . . . . . . . .
3.1 Derivation of the Analytic Gaussian equations . . . . . .
3.1.1 N-electron pulse in the absence of external forces
3.1.2 Lenses and RF cavities . . . . . . . . . . . . . . .
3.1.3 Photoemission gun . . . . . . . . . . . . . . . . .
3.1.4 Relativistic transformations . . . . . . . . . . . .
3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Drift . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Photoemission gun . . . . . . . . . . . . . . . . .
3.3.3 Lens . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.4 RF cavity . . . . . . . . . . . . . . . . . . . . . .
3.3.5 Column . . . . . . . . . . . . . . . . . . . . . . .

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

46
47
47
56
60
62
70
71
72
76
77
82
82

vi

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

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

.
.
.
.
.

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

.
.
.
.
.

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

.
.
.
.
.

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

.
.
.
.
.

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

.
.
.
.
.

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

.
.
.
.
.

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

.
.
.
.
.

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

.
.
.
.
.

.
.
.
.
.

ix

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

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

Chapter 4 Microscopic origin of the properties of TaS2
4.1 Physics of correlated electron systems . . . . . . . . . .
4.2 Structural and electronic properties of TaS2 . . . . . .
4.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Ground state of TaS2 . . . . . . . . . . . . . . . . . . .
4.4.1 Band structure and DOS . . . . . . . . . . . . .
4.4.2 Role of spin-orbit coupling . . . . . . . . . . . .
4.4.3 Role of Hubbard repulsion . . . . . . . . . . . .
4.4.4 Role of stacking disorder . . . . . . . . . . . . .
4.5 Phonons . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6 Optical absorption . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

88
88
91
95
98
98
101
103
104
108
111

Chapter 5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
Appendix A: Generation of electron coordinates . . . . . . . . . . . . . . . . . . . 122
Appendix B: Overview of Density Functional Theory . . . . . . . . . . . . . . . . 128
BIBLIOGRAPHY

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

vii

LIST OF TABLES

Table 2.1:

Table 3.1:

Table 3.2:

Table 4.1:

Table 4.2:

Physical parameters used in the simulation of the photoemission process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

Parameters used in the simulation of the UEM system using the Analytic Gaussian model. . . . . . . . . . . . . . . . . . . . . . . . . .

84

Focal spot position and corresponding pulse width for varying number
of electrons Ne and transverse photoexciting laser profile. . . . . . .

87

Theoretical and experimental (from [1]) values for the parameters of
the C-CDW unit cell . . . . . . . . . . . . . . . . . . . . . . . . . .

99

Ground state energy E0 and density of states at the Fermi energy
n(EF ) in the C-CDW and metallic phases. . . . . . . . . . . . . . . 101

viii

LIST OF FIGURES

Figure 1.1:

Figure 1.2:

Figure 2.1:

Figure 2.2:

Figure 2.3:

Spatial and temporal timescales of the various processes of interest
in Ultrafast Science. The plot was compiled with the data reported
in [2–25] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

Schematic diagram of the Ultrafast Electron Microscope system. (Image courtesy of the Ruan Group, unpublished) . . . . . . . . . . . .

5

(a) Schematic of the three-step photoemission model: (1) the electron
absorbs a photon of energy ω, (2) migrates to the surface and (3)
escapes to the vacuum. (b) Definition of the coordinate system used
in the simulation. Also shown is the Gaussian profile of the laser pulse
on the photocathode surface. The electron is emitted with velocities
given by equations (2.11)-(2.13). (c) Diagram showing the velocity of
the electron inside the metal and in vacuum. Notice the conservation
of transverse velocity across the interface and the decrease of the
velocity in the perpendicular direction, vz , due to the potential barrier
at the surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

Effect of varying the number of electrons in each macroparticle for
Ne0 = 106 , extraction field Fa =0.32 MV/m and Nmacr = 1, 10, 100.
The pulse center of mass position (¯
x, z¯), width (σx , σz ), momentum
standard deviation (σpx , σpz ) and emittance ( n,x , n,z ) are plotted
as a function of time in both the transverse (x) and longitudinal (z)
directions. Because of the radial symmetry, the x and y directions
are equal so only the parameters in the x direction are plotted. . . .

16

Time evolution of key pulse parameters, for Ne0 = 7 · 106 , extraction
fields of Fa =0.32 MV/m (continuous lines) and Fa =10 MV/m (dashed
lines): (a) number of electrons Ne (note the logarithmic scale of both
axes). (b) Pulse center of mass distance from the surface and corresponding free particle time evolution [z = 1/2(eFa /m)t2 , dotted
lines]. (c) Transverse (black) and longitudinal (red) pulse size. (d)
Transverse (black) and longitudinal (red) normalized emittance. . .

24

ix

Figure 2.4:

Figure 2.5:

Figure 2.6:

Figure 2.7:

Figure 2.8:

(a)-(d) Color map of the charge distribution in the rest frame of the
bunch projected onto the x-z plane, overlaid with arrows representing
the average electron velocity. Panel (a) shows the initial stage of
photoemission, which is independent of extraction field. Panel (b),
(c) show the intermediate and final stages for weak (Fa =0.32 MV/m)
extraction field at times t=60 ps and t=120 ps. Panel (d) corresponds
to the final stage for the strong field case (Fa =10 MV/m). Note
that the number of emitted electrons is the same in all cases, Ne0 =
7 · 106 , while the final number of emitted electrons depends on the
extraction field. (e) Transverse charge density profiles under different
extraction fields (Fa =0.32, 1 and 10 MV/m) fitted with a combination
of Gaussian (G) and Ellipsoidal (E) functions [26]. . . . . . . . . .

26

Comparison to experimental data. Panel (a) shows the raw experimental images, while panels (b)-(d) show a comparison of the longitudinal pulse profile at various times in experiment (red circles) and
simulation (black line), corresponding to an extraction field of Fa =
0.32 MV/m and Ne0 = 108 [26]. . . . . . . . . . . . . . . . . . . . . .

28

Number of electrons emitted at t=120 ps as a function of number of
generated electrons for various extraction fields. The lines serve as
a guide to the eye for the scaling Neemit ∝ Ne0 and Neemit ∝ (Ne0 )1/3 .
The data shows evidence of virtual cathode formation. Inset: Critical
number of electrons for virtual cathode formation as a function of
extraction field Fa . The dashed line shows the critical number of
electrons as a function of Fa obtained by using Eq. 2.35 [26]. . . . .

30

Longitudinal pulse width, σz , as a function of number of electrons
at time t=120 ps, for various extraction fields Fa from the simulation
(left panel) and compared to experimental data (right panel). The
lines serve as a guide to the eye for the cases σz ∝ Neemit and σz ∝
(Neemit )1/2 [26]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

(a),(b) Color map of the charge distribution in the rest frame of the
bunch projected onto the x-z plane, overlaid with arrows representing
the average electron velocity. Panel (a) shows the distribution at
t=120 ps for extraction field Fa = 10 MV/m and number of emitted
electrons Ne0 = 7 · 106 . Panel (b) shows the charge distribution under
the same conditions but with the image charge pinned to a small area
on the photocathode surface. (c) Transverse charge density profiles
under different extraction fields (Fa =0.32, 1 and 10 MV/m) showing
the central focusing peak caused by the pinned image charge. . . . .

32

x

Figure 2.9:

(a) Number of electrons emitted at t=120 ps as a function of number
of generated electrons for various extraction fields in the presence of
a pinned image charge. The red lines show the data corresponding
to the expanding image charge. (b) Longitudinal pulse width, σz ,
as a function of number of electrons at time t=120 ps, for various
extraction fields Fa with a pinned image charge. The red lines, again,
show the data corresponding to the expanding image charge. . . . .

33

Figure 2.10: (a)-(d) Color map of the charge distribution in the rest frame of the
bunch projected onto the x-z plane, overlaid with arrows representing
the average electron velocity for the top-hat (uniform) pulse. Panels
(a) and (b) show the initial and final stages of photoemission in the
presence of an image charge on the surface. Panels (c), (d) correspond
to the case without image charge. Note that the number of emitted
electrons is the same in all cases, Ne0 = 1·107 , with an extraction field
Fa = 1 MV/m . (e) Transverse charge density profiles under different
extraction fields (Fa =0.32, 1 and 10 MV/m), showing the effect of the
image charge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

Figure 2.11: (a)-(d) Color map of the charge distribution in the rest frame of the
bunch projected onto the x-z plane, overlaid with arrows representing
the average electron velocity for the elliptically shaped pulse. Panels
(a) and (b) show the initial and final stages of photoemission in the
presence of an image charge on the surface. Panels (c), (d) correspond
to the case without image charge. Note that the number of emitted
electrons is the same in all cases, Ne0 = 1·107 , with an extraction field
Fa = 1 MV/m. (e) Transverse charge density profiles under different
extraction fields (Fa =0.32, 1 and 10 MV/m), showing the effect of the
image charge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

Figure 2.12: (a) Transverse normalized emittance x and (b) longitudinal normalized emittance z as a function of number of emitted electrons Neemit
and extraction field Fa at 120 ps for a Gaussian pulse. Insets: time
dependence of x and z for three selected cases at Neemit = 107 [26].

38

Figure 2.13: (a) Transverse normalized rms emittance x and (b) longitudinal normalized rms emittance z dependence on number of emitted electrons
Neemit and extraction field Fa for varying shapes of the exciting laser
pulse at t = 120ps. Inset: number of emitted electrons Neemit at
t=120 ps as a function of Ne0 , the initial number of electrons for Fa
= 1 MV/m and varying laser profile. . . . . . . . . . . . . . . . . . .

40

xi

Figure 2.14: Time dependence of the 6D normalized emittance for Fa = 1 MV/m
both (a) below and (b) above the virtual cathode limit. The data
shown corresponds to the Gaussian (G), Elliptical (E) and Uniform
(U) cases, with the subscript 0 indicating the absence of the image
charge (G0, E0, U0 respectively). . . . . . . . . . . . . . . . . . . .

41

Figure 2.15: Time dependence of the 6D normalized emittance for Fa = 1 MV/m
both (a) below and (b) above the virtual cathode limit. The data
shown corresponds to the Gaussian (G) and Elliptical (E) cases,
for different ways of considering the image charge field: expanding,
pinned to the surface or mirroring the electron pulse distance from
the surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

Figure 2.16: (a) 4D brightness and (b) 6D brightness as a function of number
of electrons emitted, Neemit for the Gaussian, Uniform and Elliptical
cases. The data shown corresponds to Fa = 1 MV/m (closed circles)
and Fa = 10 MV/m (open circles). . . . . . . . . . . . . . . . . . . .

44

Figure 2.17: Coherence length (a) and energy spread (b) as a function of number
of electrons emitted, Neemit for the Gaussian, Uniform and Elliptical
profiles. The data shown corresponds to Fa = 1 MV/m (closed circles)
and Fa = 10 MV/m (open circles). . . . . . . . . . . . . . . . . . . .

45

Figure 3.1:

Figure 3.2:

Figure 3.3:

Figure 3.4:

Schematic representation of the pulse parameters used in the Analytic
Gaussian model in phase space. . . . . . . . . . . . . . . . . . . . .

50

Functions L(ξ) (Eq. 3.39), LT (ξ) (Eq. 3.37) and Lz (ξ) (Eq. 3.38) as
a function of ξ = σz /σT . . . . . . . . . . . . . . . . . . . . . . . .

54

Effect of stepwise electric field of photoemission gun: (a) z < d, (b)
z >d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

Color map of the charge distribution of an electron bunch below the
virtual cathode limit (Fa =1MV/m, Ne =106 , Gaussian laser profile)
in the rest frame of the bunch projected onto the px -x plane [(a),
(b)] and pz -z plane [(c),(d)]. Also shown is the time evolution of the
Analytic Gaussian model for the pulse, with the black line indicating
the contour where the intensity drops to 0.5 times its peak value. . .

73

xii

Figure 3.5:

Time evolution of pulse parameters below the virtual cathode limit
(Fa =1MV/m, Ne =106 , Gaussian laser profile) comparing N-particle
simulations and the Analytic Gaussian model. Panel (a) shows the
pulse center of mass position, panels (b)-(d) show the transverse (T)
√
and longitudinal (z) pulse width σi [panel (b)], chirp γi [panel (c)]
√
and momentum spread σpi [panel (d)]. . . . . . . . . . . . . . . . .

74

Color map of the charge distribution of an electron bunch below the
virtual cathode limit (Fa =1MV/m, Ne =106 , Elliptical laser profile)
in the rest frame of the bunch projected onto the px -x plane [(a),
(b)] and pz -z plane [(c),(d)]. Also shown is the time evolution of the
Analytic Gaussian model for the pulse, with the black line indicating
the contour where the intensity drops to 0.5 times its peak value. . .

75

Time evolution of pulse parameters below the virtual cathode limit
(Fa =1MV/m, Ne =106 , Elliptical laser profile) comparing N-particle
simulations and the Analytic Gaussian model. Panel (a) shows the
pulse center of mass position, panels (b)-(d) show the transverse (T)
√
and longitudinal (z) pulse width σi [panel (b)], chirp γi [panel (c)]
√
and momentum spread σpi [panel (d)]. . . . . . . . . . . . . . . . .

76

Color map of the charge distribution of an electron bunch above the
virtual cathode limit (Fa =1MV/m, Ne =108 , Gaussian laser profile)
in the rest frame of the bunch projected onto the px -x plane [(a),
(b)] and pz -z plane [(c),(d)]. Also shown is the time evolution of the
Analytic Gaussian model for the pulse, with the black line indicating
the contour where the intensity drops to 0.5 times its peak value. . .

77

Time evolution of pulse parameters above the virtual cathode limit
(Fa =1MV/m, Ne =108 , Gaussian laser profile) comparing N-particle
simulations and the Analytic Gaussian model. Panel (a) shows the
pulse center of mass position, panels (b)-(d) show the transverse (T)
√
and longitudinal (z) pulse width σi [panel (b)], chirp γi [panel (c)]
√
and momentum spread σpi [panel (d)]. . . . . . . . . . . . . . . . .

78

Figure 3.10: Longitudinal pulse parameters as a function of distance along the
microscope column for Fa =1MV/m, as the pulse crosses the anode
plane (at z=45mm, indicated with dashed line). Panel (a) shows the
center of mass velocity as a fraction of the speed of light c, while
panels (b)-(d) show the variation of the longitudinal (z) pulse width
√
σz [panel (b)], chirp γz /γ0 , where γ0 is the initial chirp [panel (c)]
√
and local momentum variance ηz [panel (d)]. . . . . . . . . . . . .

79

Figure 3.6:

Figure 3.7:

Figure 3.8:

Figure 3.9:

xiii

Figure 3.11: (a) Fractional lens strength M (z)/M0 as a function of position. (b)(d) Transverse pulse parameters as a function of distance along the
microscope column for Ne =106 , with a magnetic lens of varying strength
M0 at z=45mm (indicated with a dashed line). Shown are the trans√
verse (T) pulse width σT [panel (b)], chirp γz /γ0 [panel (c)] and
√
local momentum variance ηz [panel (d)]. . . . . . . . . . . . . . . . 80
√
Figure 3.12: Transverse (T) pulse width σT as a function of distance along the
microscope column with a magnetic lens (M0 = 5) at z=45mm (indicated with a dashed line) for varying numbers of electrons in the
pulse. (a) Effect of lens given the same initial conditions. (b) Effect
of lens using the initial conditions from the photoemission simulations. 81
Figure 3.13: Pulse parameters as a function of distance along the microscope column for Ne =106 with an RF cavity centered at z=45mm (indicated
by the vertical dashed line) for different RF field strengths. (a) Center
of mass velocity as a fraction of speed of light c and electric field as a
fraction of the peak value E0 . (b)-(d) Transverse (T) and longitudinal
√
(z) pulse width σi [panel (b)], chirp γi /γ0 [panel (c)] and local mo√
mentum variance ηi [panel (d)]. The continuous (dashed) lines indicate the longitudinal (transverse) components and the colors correspond, respectively, to Erf = 2.4M V /m (black) and Erf = 1.0M V /m
(red). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83

√
√
Figure 3.14: (a) Transverse pulse width σT and (b) longitudinal pulse width σz
as a function of distance along the microscope column for varying Ne ,
with an RF cavity centered at z=45mm (indicated by the vertical
dashed line) with Erf =2.4MV/m. . . . . . . . . . . . . . . . . . . .

84

Figure 3.15: Longitudinal (black) and transverse (dashed red) pulse parameters
as a function of distance along the microscope column for Ne =106 .
√
Shown are the (a) pulse width σi , (b) chirp γi and (c) local momen√
tum variance ηi . The region corresponding to the photoemission
gun is shown in blue, while that corresponding to the RF cavity is in
green. The dashed lines correspond to the magnetic lenses. . . . . .

85

Figure 3.16: Detail of the transverse and longitudinal pulse widths in the optimized
column around the focal spot for a Gaussian laser pulse with (a)
Ne = 105 , (b) Ne = 106 and (c) Ne = 107 electrons and (d) for an
elliptical laser pulse with Ne = 106 electrons. . . . . . . . . . . . . .

86

xiv

Figure 4.1:

Figure 4.2:

Figure 4.3:

Figure 4.4:

Figure 4.5:

Figure 4.6:

Figure 4.7:

Figure 4.8:

Three classes of insulators and their respective characteristic response
times. (a) Mott insulator (b) Excitonic insulator (c) Peierls insulator.
In each case the gray shading and dashed lines indicate the metallic state and the solid lines and colored shading show the resulting
insulating state both in real and momentum space. (d) Timescales
corresponding to probing these processes with a pump-probe experiment. [27] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90

TaS2 structure corresponding to the primitive unit cell (undistorted
metallic phase): (a) 3D view, (b) Top view and (c) Side view. Ta
atoms are drawn in blue, S atoms in red. The black lines delimit the
primitive 1x1 unit cell. . . . . . . . . . . . . . . . . . . . . . . . . .

92

Positions of the Ta atoms corresponding to the reconstructed “Star of
David” arrangement in a 2x2 unit cell. The undistorted positions are
shown in yellow, while the reconstructed lattice and bonds resulting
from the contraction associated with CDW formation are drawn in
blue. The labels correspond to the three inequivalent atom types of
the distorted cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93

Resistivity as a function of temperature in 1T-TaS2 . The phases
shown are: C-CDW below 190K, NC-CDW below 350K, IC-CDW
below 550K and metallic for T>550K. Also shown is the trigonal
phase present during the heating cycle between 200-280K. [28] . . .

93

Convergence of the ground state energy for the metallic phase as a
function of energy cutoff and k point mesh size. The dashed lines
correspond to the acceptable range of values that gives a precision of
1meV/atom. Note that due to the size of the reciprocal unit cell, the
Nth k-point mesh corresponds to a mesh of size N in the kx and ky
directions and N/2 in the kz direction. . . . . . . . . . . . . . . . . .

96

Convergence of the ground state energy for the distorted phase as a
function of energy cutoff and k point mesh size for varying energy
cutoff. The dashed lines correspond to the acceptable range of values
that gives a precision of 1meV/atom. Note that due to the size of
the reciprocal unit cell, the Nth k-point mesh corresponds to a mesh
of size N in the kx and ky directions and 2 · N in the kz direction. .

96

Top (a) and side (b) views of the relaxed C-CDW structure, with
the black line indicating the boundary of the unit cell. Ta atoms are
drawn in blue, while S atoms are in red. . . . . . . . . . . . . . . . .

99

First Brillouin zone corresponding to the C-CDW TaS2 unit cell. The
labels represent the conventional high symmetry points. [29] . . . . .

99

xv

Figure 4.9:

Band structure along high symmetry lines and density of states for
the C-CDW phase. The contribution to the DOS due to the Ta d
orbitals and that due to only the dz2 is also plotted. . . . . . . . . . 100

Figure 4.10: Band structure along high symmetry lines and density of states for
the metallic phase. The contribution to the DOS due to the Ta d
orbitals and that due to only the dz2 is also plotted. . . . . . . . . . 101
Figure 4.11: Band structure along high symmetry lines and density of states for
the C-CDW phase both with (black) and without (red) spin orbit
coupling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
Figure 4.12: Detail of the density of states for the C-CDW phase both with (continuous lines) and without (dashed lines) spin orbit coupling. (a)
Atom projected densities. (b) Orbital projected densities. . . . . . . 103
Figure 4.13: Band structure along high symmetry lines and density of states for
the metallic phase both with (black) and without (red) spin orbit
coupling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
Figure 4.14: Band structure along high symmetry lines and density of states for
the C-CDW phase in the DFT+U approximation with varying values
of the on-site electron-electron interaction U . . . . . . . . . . . . . 105
Figure 4.15: Units cells for the stacking sequences considered: (a) A stacking, (b)
C stacking and (c) AC stacking. . . . . . . . . . . . . . . . . . . . . 106
Figure 4.16: Band structure along high symmetry lines and density of states for the
C-CDW phase with different stacking: (a) A stacking, (b) C stacking
and (c) AC stacking. The last two panels also show a comparison to
the DOS obtained with the A stacking (stacking without disorder). . 107
Figure 4.17: Calculated phonon band structure for the metallic (a) and C-CDW
(b) phases. Imaginary frequencies, corresponding to unstable modes
are represented by negative values. (c) Phonon density of states for
the metallic (green) and C-CDW (black) phases. The continuous line
represents the total DOS, while the dashed line shows the contribution due to the Tantalum atom vibrations. . . . . . . . . . . . . . . 109
Figure 4.18: Detail of the phonon density of states per TaS2 for the C-CDW phase
in the range 0 meV to 18 meV (corresponding to 0 THz to 4 THz). The
red lines indicate the site projected DOS. . . . . . . . . . . . . . . . 111

xvi

Figure 4.19: In-plane optical parameters for TaS2 : real ( 1 ) and imaginary ( 2 )
parts of the dielectric constant for the metallic (a) and C-CDW (b)
phases. The insets show a detail of the dielectric constant in the
energy range relevant for pump-probe experiments. (c) Absorption
coefficient (d) Energy-loss function . . . . . . . . . . . . . . . . . . . 114
Figure 4.20: Out of plane optical parameters for TaS2 : real ( 1 ) and imaginary
( 2 ) parts of the dielectric constant for the metallic (a) and C-CDW
(b) phases. The insets show a detail of the dielectric constant in the
energy range relevant for pump-probe experiments. (c) Absorption
coefficient (d) Energy-loss function . . . . . . . . . . . . . . . . . . . 115

xvii

Chapter 1
Introduction

1.1

Why study ultrafast science

Phenomena in the natural world take place on a variety of time scales, ranging from the
millions of years one deals with in cosmology, to the infinitesimally small fractions of seconds
involved with chemical, atomic and nuclear reactions. As the title of this thesis suggests,
in the present work we are interested in ultrafast science, that is we seek to understand the
transient non-equilibrium processes that happen in materials on timescales of less than 10−9
seconds. Very broadly speaking, there is also often a correlation between timescale and size in
physical systems, meaning that as the time scales of the phenomena we are interested in get
smaller, so do the length scales over which these take place. For this reason when speaking
of ultrafast science, we are also implicitly talking about length scales roughly between 10−6
and 10−10 meters.
With these considerations in mind, we have a very broad field that includes a number
of interesting time dependent processes ranging from biology [2–5] to material science (from
ultrafast phase transitions [6–11] to charge transport at nanointerfaces [12,13] and dynamics
1

in correlated electron materials [14–16]), from chemistry [17–20] to plasma physics [21–24].
As seen in the examples of Fig. 1.1 these phenomena take place across several different
orders of magnitude of spatial and temporal resolution, but in each case further understanding of the fundamental interactions that determine the system’s properties requires gaining
more insight into atomic motion on a timescale corresponding to the period of the atomic
vibrations. In addition to this, further understanding in this field would open the doors for
numerous technological applications [30].

Figure 1.1: Spatial and temporal timescales of the various processes of interest in Ultrafast
Science. The plot was compiled with the data reported in [2–25]

1.2

Pump probe microscopy

In practice these systems are best explored by using pump-probe measurements: we excite
the system under consideration by perturbing it with a pump and measure its response at
2

various delay times with a probe. By recording the response of the system as a function
of delay time, we can obtain a description of the time evolution of the sample. In order to
resolve these ultrafast phenomena, both the pump and the probe pulse need to be sufficiently
short in time, which is the reason why the vast majority of this type of measurement use
a femtosecond laser [31–33] to excite the system in virtue of its short duration, high peak
power, broad bandwidth and commercial availability.
As far as probes are concerned, the main candidates are x-rays and electrons [25, 34, 35].
The advantage of using x-rays is that they offer Angstrom size wavelength, high spatial
coherence and short pulse length, but their generation with a high enough brightness requires
large scale facilities (such as SLAC’s Linac). In addition to this, the interaction of x-rays with
matter has a very small cross section, meaning that a high intensity and sensitive detection
are required to achieve a good signal to noise ratio.
Electron probes, on the other hand, are easily generated using tabletop scale equipment,
deal less damage to the sample compared to x-rays and have a high scattering cross section.
As a consequence of this high cross section, the mean free path of the probe electrons inside
the material is shorter than for x-rays and typically better matches the optical penetration
depth, allowing to only probe the pumped volume. The drawback presented by this probe is
that one has to deal with the space charge broadening due to the Coulomb repulsion between
electrons, which leads to a decreased coherence length and increased pulse size. In addition,
due to the short penetration length, there is a limited thickness that can be measured making
experiments in liquid phases very challenging.
Given these considerations, one should view electrons and x-rays as complementary
probes to use in ultrafast studies of materials. For the present work, though, we will only
discuss systems using electrons as a probe.
3

1.3

Ultrafast Electron Microscope at MSU

In the Ultrafast Research group at Michigan State University, we are working on the development of an Ultrafast Electron Microscope (UEM) that would allow measurements of some
of the processes described in the previous sections. More specifically, we want to perform
single-shot diffraction to observe irreversible processes in materials, such as photo-induced
non equilibrium dynamics of lattice and charge transfer. Single-shot means that enough
electrons for good signal to noise ratio are contained in a single probe pulse, in contrast to
the stroboscopic mode in which the measurement for a given time delay between pump and
probe is repeated over many shots to collect enough signal. The Rose criterion tells us that
approximately 100 electrons per pixel are needed to resolve a pattern, which translates to
roughly 105 to 107 electrons needed in a single pulse for diffraction and 107 to 109 electrons
if one is interested in imaging.
The requirement of such a high number of electrons means that the expansion of the
pulse due to space charge effects must be carefully considered and reversed in order to
focus the pulse at the sample. In the transverse direction this can be easily done by using
magnetic lenses but focusing in the longitudinal direction is harder to achieve. The solution
we adopted is to use a radio frequency (RF) cavity phase locked with the photoexciting
laser that generates the electron pulse. A schematic of the UEM system is shown in Fig.
1.2 and further details can be found in recent dissertations by Zhensheng Tao and Kiseok
Chang. which cover the experimental development. This thesis covers some of the theoretical
components of the project.

4

Figure 1.2: Schematic diagram of the Ultrafast Electron Microscope system. (Image courtesy
of the Ruan Group, unpublished)

1.4

Overview

The central focus of this dissertation is to gain more insight into the development of the UEM
system and on its applications. To do so, in the first part of this work, we give an overview
of the simulations relative to the UEM system. In Chapter 2, we look at the details of the
5

photoemission process under different experimental conditions and discuss how to generate
optimal electron pulses for our application. The code used for the simulations was developed
by He Zhang [36]; using it to understand the photoemission process and extending it to
different laser and surface image charge configurations was my own contribution.
In Chapter 3 we proceed to simulations of the whole microscope column using a mean
field model and provide estimates of the achievable resolution in the system, given different
initial conditions. This so-called Analytic Gaussian model was originally derived by Michalik
and Sipe [37, 38] and refined by Berger and Schroeder [39]; my contribution was extending
it to treat the photoemission gun and the relativistic motion of the electrons and writing a
program to apply it to the study of the UEM system and combine it with an optimization
routine.
Chapter 4 moves away from microscope design and presents an example of a system,
TaS2 , displaying interesting ultrafast phenomena that can be probed using the microscope.
By performing ab-initio calculations of the band structure, phonon modes and optical properties of this material, we shine light on the microscopic mechanisms underlying the insulating
gap seen in experiments. All the simulations in this chapter were performed by myself using
the well established Vienna Ab-initio Simulation Package (VASP).
Chapter 5 presents the conclusions of this thesis and an outlook on future research
directions.

6

Chapter 2
Generation of the electron pulse

2.1

Importance of the initial stages of pulse generation

Once an electron pulse has been generated in the initial stage of the UEM microscope
introduced in the previous chapter, it can be shaped as it travels through the column using
a variety of technologies such as magnetic lenses and radio frequency cavities [5, 25, 40–45].
What sets limits on the achievable temporal and spatial resolutions, though, is the phase
space volume the pulse occupies, which in turn is defined by the photoemission conditions
and the initial Coulomb expansion.
In the transverse direction our ability to resolve spatial features can be quantified with
the coherence length Lc , which at the beam waist can be expressed as [46]

Lc =

h σx
,
2m0 c n,x

(2.1)

where we have assumed that the transverse and longitudinal directions are decoupled, h is
Plank’s constant, m0 is the rest mass of the electron, c is the speed of light, σx is the beam

7

radius and

n,x

is the normalized transverse emittance defined as [47–53]

n,x

=

1
mc

x2 p2x − xpx 2 .

(2.2)

In the longitudinal direction, using the uncertainty principle and assuming a Gaussian
pulse shape, we have that
∆E =
where

n,z

z

∆t

γm0 c.

(2.3)

is the normalized longitudinal emittance, ∆t and ∆E are, respectively, the tem-

poral and energy resolutions and γ is the relativistic Lorentz factor.
One last figure of merit for electron sources is the normalized beam degeneracy, defined
as [40, 43]
η = B6D
where B6D is the 6D beam brightness and

0

3
0

=

Ne
x y z

3
0,

(2.4)

= h/(m0 c) is the emittance quantum. Consider-

ing that a single electron can have two possible spin orientations, the upper limit calculated
using the Pauli exclusion principle is η = 2, with typical values ranging from
cold field emission gun to

10−4 for a

10−12 if using a thermionic gun [40]. An ideal source of electron

pulses would maximize the coherence length Lc and the degeneracy η, while maintaining a
low energy spread ∆E which, in turn, requires minimizing the longitudinal emittance

n,z .

It is worth noting here that the emittance defined in Eq. 2.2 is a statistical approximation
of the true 6D emittance, where correlations between the x, y and z directions have been
neglected. We calculate this quantity and parameters derived from it as traditionally the
emittance has been used as a key figure of merit for UEM systems, being an estimate
of the phase space occupied by the pulse and therefore setting limits on the achievable

8

resolution [26].
For all these reasons it is very important for the design of high brightness, high resolution
electron beams to control the initial stages of pulse formation, so that the increase in beam
emittance due to space charge effects can be minimized while maintaining a high number
of electrons. To better understand the nonlinear interplay between the space charge driven
expansion and laser and photocathode material properties on the generated electron bunch,
we have conducted explicit N-particle simulations of the electron photoemission process
using COSY Infinity, a code designed for high performance scientific computing and beam
dynamics simulations [54]. In the following sections we will first introduce the model used
to describe this process, followed by a description of the key elements of the simulation and
then proceed to a discussion of the results and their comparison with experimental data.

2.2

Modeling the photoemission process

In the UEM system, the electron bunch is generated through photoemission from a gold
photocathode irradiated with a 50fs laser pulse. This process can be described using the socalled three-step photoemission model [48,55] in which each electron is emitted independently
as a result of (1) absorbing a photon of energy ω, (2) diffusing to the surface and (3) escaping
to the vacuum [Fig. 2.1(a)]. As the thickness of the photocathode used in the experiment
is smaller than the electron mean free path, electron-electron scattering does not contribute
and step (2) can be neglected by simply treating all of the electrons as generated on the
surface. We will now look at the remaining steps in the photoemission process in more
detail.
The photocathode is irradiated with a laser pulse with energy ω=4.66 eV. As the photon

9

energy is higher than the gold work function, W =4.0 eV to 4.6 eV [41, 56–58], a fraction of
the electrons inside the metal can absorb enough energy to overcome the potential barrier
and escape to the vacuum. The energy difference, ∆E = ω − W , allows only electrons with
energies in the range [EF −∆E, EF ] to be emitted, and since ∆E is small, the generated bunch
has a narrow energy spread. The lower limit for this range is set by the condition that the
energy of the electron after the collision with the photon must satisfy E = Ei + ω > EF +W ,
while the upper limit reflects the fact that only states up to the Fermi energy are initially
occupied. We assume an uniform distribution of states within this energy band.
In order to escape the photocathode, the electrons need to cross the metal-vacuum boundary with sufficient velocity in the direction perpendicular to the surface:

2(EF + W )
,
m

vz ≥ vz min =

(2.5)

where vz min is the velocity corresponding to the minimum kinetic energy required to overcome the potential barrier EF + W and the z axis is assigned to the direction perpendicular
to the metal surface [Fig. 2.1(b)].
Given an electron with initial energy Ei ∈ [EF − ∆E, EF ], we have that

vi =

2(Ei + ω)
,
m

(2.6)

which, together with the equation for vz min (Eq. 2.5), determines a maximum internal angle
for which the electron can escape:

cos θi,max =

vz min
=
vi

10

EF + W
.
Ei + ω

(2.7)

E
EF +W
EF
F- E

metal

vacuum

z

Figure 2.1: (a) Schematic of the three-step photoemission model: (1) the electron absorbs
a photon of energy ω, (2) migrates to the surface and (3) escapes to the vacuum. (b)
Definition of the coordinate system used in the simulation. Also shown is the Gaussian
profile of the laser pulse on the photocathode surface. The electron is emitted with velocities
given by equations (2.11)-(2.13). (c) Diagram showing the velocity of the electron inside the
metal and in vacuum. Notice the conservation of transverse velocity across the interface and
the decrease of the velocity in the perpendicular direction, vz , due to the potential barrier
at the surface.
Therefore inside the cathode:

vx in = vi sin θi cos φi ,

(2.8)

vy in = vi sin θi sin φi ,

(2.9)

vz in = vi cos θi ,

(2.10)

where θi ∈ [0, θmax ] and φi ∈ [0, 2π]. Outside the cathode the velocities in the directions
parallel to the surface are conserved, while the velocity in the z direction is reduced due to

11

the potential barrier [see Fig. 2.1(c)]:

vx out = vx in = vi sin θi cos φi ,

(2.11)

vy out = vy in = vi sin θi sin φi ,

(2.12)

vz out =

vz2 in − 2(EF + W )/m =

(vi cos θi )2 − 2(EF + W )/m.

(2.13)

In addition to defining the initial velocities of the electrons with the above equations,
the initial positions must also be calculated. The coordinates in directions parallel to the
surface (x, y in our notation) are randomly generated from the spatial profile of the laser
pulse. More details on how this is done in practice for the pulse profiles considered are given
in Appendix 5. The initial z coordinate is set equal to zero since, as previously mentioned,
we treat all of the electrons as generated on the surface. The physical parameters used to
describe the system in the simulation are summarized in Table 2.1.
Fermi energy, EF
5 eV
Work function, W
4.45 eV
Photon energy, ω
4.66 eV
Laser pulse duration, ∆t
50 fs
Laser pulse width, σr
91 ➭m
Table 2.1: Physical parameters used in the simulation of the photoemission process.

In each simulation run, given the total number of electrons to generate, Ne0 , the number of
electrons emitted during each infinitesimal time interval can be calculated assuming that the
√
laser pulse has a Gaussian distribution in time with a standard deviation σ = ∆t/(2 2 ln 2),
where ∆t is the laser duration (FWHM). Considering then the range of times [−3σ, +3σ]
and dividing it into N parts (typically N = 100), it is possible to calculate the number of

12

electrons generated in each ith interval as
√
√
2 ) − erf(t
erf(t
/
2σ
/
2σ 2 )
i
i−1
,
Ni = Ne0
2C

(2.14)

where erf is the error function and

1
C= √
σ 2π

+3σ

x2

dxe− 2σ2

0.997.

(2.15)

−3σ

Note that the quantity Ne0 indicates the number of electrons that are emitted due to the laser
illumination and is therefore directly proportional to the laser fluence F . The proportionality
constant depends on the quantum efficiency as:

F =

ωNe0
QE · (1 − R),
A

(2.16)

where A = πσx σy is the area of the laser spot on the cathode surface, QE is the quantum
efficiency, typically on the order of

10−4 [48], and R is the reflectivity of the sample surface.

To conclude the discussion on the model used, it is necessary to point out that each
emitted electron contributes to the formation of an image charge on the surface of the
photocathode model. This positive field plays an important role in shaping the electron
pulse dynamics and will be discussed in more detail in section 2.4. Taking this into account,
each electron in the system is treated as subject to both the space charge field between
the particles in the bunch and the positive field generated by the image charge on the
surface. In addition to the Coulomb electron-electron forces and the field due to the image
charge, a constant electric field is applied in the direction perpendicular to the photocathode
surface. Typical values for extraction fields are between 0.1-1MV/m for the thermionic gun
13

geometry [40], between 1-10MV/m for the DC electron gun [59] and 10-100MV/m for the
RF gun [60]. To provide relevant data over the whole range, we vary the extraction field
between 0.1MV/m and 10MV/m which corresponds to generation of electrons with energies
ranging from 5 · 10−3 keV to 100 keV at t=120ps.

2.3

Details of the simulations

As we wish to simulate the time evolution of the electron bunch after the photoemission
process has taken place, it is necessary to find a good scheme to treat the Coulomb interaction
that would be both accurate and computationally advantageous. In the present work we have
used COSY INFINITY [54], which treats the space charge effects with a Multiple Level Fast
Multipole Method (MLFMM). In short, the MLFMM uses a tree algorithm to partition the
simulated space into boxes until the number of charges in each one is lower than a certain
threshold (typically NBOX < 100), after which the near and far region of each box are
determined. The near region includes the adjacent boxes, for which the electromagnetic field
is calculated directly. For the far region, a multipole expansion is used to calculate the field
in the box of interest. This approach gives an algorithm that scales almost linearly with
the number of particles with a prefactor that depends on the multipole order used. Further
detail on the MLFMM are given elsewhere [36, 61].
In order to further reduce the computational time and allow simulations with higher numbers of electrons, we use macroparticles so that typically one simulated particle corresponds
to 100 electrons. The effect of using macroparticles made of different numbers of electrons
(NM P = 1, 10, 100) was checked and the results for a typical set of parameters are shown
in Fig. 2.2. The use of macroparticles does not affect the simulation results significantly,

14

with the biggest effect being a constant shift of about 6% of the transverse momentum standard deviation, σpx for the NM P = 100 case. This shift is propagated to the calculated
emittance

n,x

but, since it is only fractions of µm and since our main interest is comparing

values obtained with the same NM P varying the other emission conditions, it does not have
a significant impact on the discussion that follows.
To treat the time evolution of the emitted electrons, their equations of motion are solved
at each time step using a fourth order Runge Kutta algorithm. The particles that fall back
on the surface (z < 0) are removed from the simulation. Collisions between electrons are
avoided by imposing a maximum force cutoff, which is translated into a constraint on the
minimum allowed distance between two particles: rmin = γra , where γ = 0.001 and ra is the
average distance of electrons in the selected box. For an analysis of the effect of imposing
this cutoff, see Chapter 5 in [36].
The value of the time step, ∆t, is adjusted throughout the simulation in order to achieve
a good compromise between numerical accuracy and simulation length. In particular during
the initial stages when the electrons are photoemitted and are very close to each other and
to the surface, a shorter time-step is needed to keep the forces finite and avoid excessive
errors due to the force cutoff. As the electron pulse expands, progressively longer time-steps
can be used. For this reason, while the electrons are being inserted into the simulation box,
∆t = 3σt /50 = 1.26 fs. After this phase concludes (at Nstep = 100), ∆t is increased every 100
time steps to, respectively, ∆t =4.5, 14, 80fs until reaching the final value of ∆t = 0.55 ps,
which is kept constant for the remainder of the simulation.
Typical simulations run for Nstep = 620, which corresponds to a final time of t = 120 ps
and take about 1.5hrs on 8 cores, for Ne0 = 1 · 107 electrons and extraction field Fa =
1 MV/m. The standard parallel MPI routines implemented in COSY INFINITY and its
15

20
15
10
5
0
-5
-10
-15
-20

Longitudinal
400
350
300
250
200
150
100
50
0

z (µm)

Nmacr = 100
Nmacr = 10
Nmacr = 1

–

–

x (µm)

Transverse

0

20

40

60

80

100 120

Nmacr = 100
Nmacr = 10
Nmacr = 1

0

20

40

115
110
105
100
95
90

σpx (10-6 MeV/c)

500
450
400
350
300
250
200
150
100
50

20

40

60
t(ps)

σz (µm)

Nmacr = 100
Nmacr = 10
Nmacr = 1

0

80

20

40

60
t(ps)

80

40

60
t(ps)

100 120

80

100 120

10-2
10-3
10-4
10-5
10-6
10-7

20

40

60
t(ps)

80

100 120

80

100 120

Nmacr = 100
Nmacr = 10
Nmacr = 1

0

εn,z (µm)

εn,x (µm)

20

1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

100 120

Nmacr = 100
Nmacr = 10
Nmacr = 1
0

80

Nmacr = 100
Nmacr = 10
Nmacr = 1

0

10-1

10-2

70
60
50
40
30
20
10
0

100 120

Nmacr = 100
Nmacr = 10
Nmacr = 1

0

60
t(ps)

σpz (10-3 MeV/c)

σx (µm)

t(ps)

20

40

60
t(ps)

Nmacr = 100
Nmacr = 10
Nmacr = 1
0

20

40

60
t(ps)

80

100 120

Figure 2.2: Effect of varying the number of electrons in each macroparticle for Ne0 = 106 ,
extraction field Fa =0.32 MV/m and Nmacr = 1, 10, 100. The pulse center of mass position
(¯
x, z¯), width (σx , σz ), momentum standard deviation (σpx , σpz ) and emittance ( n,x , n,z )
are plotted as a function of time in both the transverse (x) and longitudinal (z) directions.
Because of the radial symmetry, the x and y directions are equal so only the parameters in
the x direction are plotted.
MLFMM extension, were used and the parallel simulations were run on the cluster of the
High Performance Computer Center (HPCC) at Michigan State University.
At each step the pulse parameters (position and momentum averages and standard de16

viations) are saved to a text file and every 10 steps the positions and momenta of all of the
electrons are saved to a separate file for post processing purposes, such as the calculation of
the normalized rms emittance (Eq. 2.2).

2.4

Treatment of the positive surface field

As briefly mentioned, the image charge generated on the surface by the emitted electrons
plays an important role in determining the final properties of the electron bunch. In order to
model this effect, the charge on the surface is approximated with a continuous distribution
and its field is calculated by dividing the surface into circular strips with a width ∆r

σr ,

where σr is the laser pulse width. Each strip is approximated as having a uniform charge
density and, thanks to the circular symmetry of our problem, the field due to each one can
be calculated using the well known equation for the field of a uniformly charged ring [62].
The total radial and longitudinal electric fields are then simply given by a sum of the fields
generated by each ring. Note that the value of the positive charge on the surface depends
on the number of electrons in the pulse and is therefore reduced by electron recombination.
The equations describing the details of this positive field depend on the shape of the exciting
laser pulse and will be presented in the next section.

2.4.1

Equations for the surface field: Gaussian, Uniform and Elliptical cases

Since the laser pulse used for the photoemission is radially symmetric in the photocathode
plane, we take a circular area with radius R0 around the center of the laser spot and divide it
into Nr = 1000 concentric rings where each has width R0 /Nr and radii ri,int = (i − 1)R0 /Nr ,
17

ri,ext = iR0 /Nr with i ∈ [1, Nr ]. The charge in each ring can now be calculated as
2π

Qi = Q

ri,ext

dr r · f (r),

dθ
0

(2.17)

ri,int

where Q = Ne e is the total charge in the system and f (r) is the normalized charge distribution on the surface.
Since R0 /Nr

σr , we now approximate each ring as being infinitesimally thin and

located at ri,0 = (ri,int + ri,ext )/2 with charge Qi . Given this, the radial and longitudinal
electric fields are simply given by the textbook equations [62] for a ring of charge as:
E(k) r2 + rr0 R
2Qi 1
−
Er (r, z) =
π rR2 1 − k 2
R
2
Ez (r, z) =

E(k)
− K(k)
1 − k2

2Qi 1 E(k)
(z − z0 ),
π R3 1 − k 2

,

(2.18)
(2.19)

where K(k) and E(k) are the complete elliptic integrals of first and second kind respectively:
π/2

dθ

K(k) =

1 − k 2 sin2 θ

0

,

(2.20)

π/2

E(k) =

dθ

1 − k 2 sin2 θ,

(2.21)

0

and

k2 =

4rr0
,
R2

(2.22)

R2 = (r + r0 )2 + (z − z0 )2 ,

(2.23)

where r0 is the radius of the ring and z0 is the coordinate of its center, with z0 = 0 in our

18

system.
The sum of the fields generated by each ring i with i ∈ [1, Nr ] is then calculated on a
mesh and saved to a file for use in the simulation. The mesh is typically a square 201 x
201 grid from 0 mm to 1.16 mm in both the longitudinal and transverse directions and the
COSY-MLFMM code linearly interpolates the electric field between the mesh points.
Details on the laser pulse shape enter in the calculation through the function f (r) that
specifies the charge distribution on the surface and is used in the calculation of the charge
of each ring Qi . For a Gaussian distribution

f (r) =

r2
1
exp
−
2πσr2
2σr2

,

(2.24)

and
2π

Qi = Q

ri,ext

dr r · f (r) = Q exp −

dθ
0

ri,int

2
ri,int
2σr2

− exp −

2
ri,ext
2σr2

.

(2.25)

We choose R0 = 4σr so that ri,int = (i−1)4σr /Nr , ri,ext = i4σr /Nr and ri,0 = 4σr /Nr (i−1/2).
In the case of a Uniform (top-hat) distribution

f (r) =





C =

1
πR02




0

r ≤ R0
(2.26)
r > R0

where R0 is the radius of the laser pulse on the surface and the constant is chosen so that
f (r)rdrdθ = 1. In this case the charge contained in each ring is given by
2π

Qi = Q

ri,ext

dr r · f (r) = Q

dθ
0

ri,int

19

ri,ext
R0

2

−

ri,int
R0

2

.

(2.27)

Since ri,int = (i − 1)R0 /Nr , ri,ext = iR0 /Nr and ri,0 = R0 /Nr (i − 1/2), we can rewrite this as

Qi = Q

2i − 1
.
Nr2

(2.28)

We choose R0 = 2σr , so that the standard deviation of the Gaussian and Uniform distributions in the transverse direction is the same.
In the elliptical case, the laser pulse has a distribution on the surface given by

3
f (r) =
2πR02

1−

r
R0

2

,

(2.29)

where R0 is again the radius of the laser pulse on the surface and 0 < r < R0 . Each ring on
the surface contains a fraction of the total charge equal to
ri,ext

2π

dr r · f (r) = Q

dθ

Qi = Q
0

ri,int

2
ri,int
1− 2
R0

3/2

2
ri,ext
− 1− 2
R0

3/2

.

Again we wish to compare distributions of equal width and therefore choose R0 =

(2.30)
√
5σr so

that the standard deviation of the transverse distribution is the same as in the Gaussian
case.

2.4.2

Time evolution of the positive field

In the treatment presented so far we implicitly assumed that the positive field on the surface
had no time dependence. In the initial stages of the simulation the shape of the electron
pulse mirrors the profile of the exciting laser pulse, so it is natural to assume that this would
be also true for the positive charge on the surface. Given R0 , the spatial width of the image
charge distribution on the surface, we can safely assume that R0 (t = 0) corresponds to the
20

width of the laser profile in the direction parallel to the surface. In general R0 does depend
on time and its temporal evolution has a complex dependence on the instantaneous shape
of the emitted pulse and on the cathode properties, as both surface characteristics and film
thickness play a role in determining the image charge created by the electrons.
To simplify this complex dependence into something treatable, the experimental systems
can be divided into roughly two categories. If the positive charges are pinned on the surface,
as would be the case for example in a patterned photocathode [63], the expansion of R0 will
be hindered or even blocked. We approximate these cases by taking R0 (t) = R0 (t = 0), so
that the image charge is pinned to a small constant area on the surface corresponding to its
initial value.
On the other hand if the positive charges on the surface are able to expand and mirror
the dynamics of the electron pulse, such as for a flat metal surface, the radius of the positive
charge R0 (t) will increase as the electron pulse spreads. This can be treated in a simple way
by taking the standard deviation of the surface charge distribution equal to the instantaneous
transverse width (standard deviation) of the electron pulse. This means that each uniform
ring of charge Qi has a radius ri,0 that increases with time.
In the simulation code the positive field is calculated on a mesh, so we note that if, for
example, we double the radius of each ring, r0 = 2r0 , and also double the coordinates of
each mesh point r = 2r, z = 2z, then we get

R 2 = (r + r0 )2 + z 2 = 4 (r + r0 )2 + z 2 = 4R2 ,
k2 =

4(2r)(2r0 )
4rr0
4r r0
=
= 2 = k2.
2
2
R
4R
R

(2.31)
(2.32)

This means that in our expression for the electric field, the terms depending on k are not

21

affected by this transformation, while the ones depending on R would be doubled. This
translates to a transformation of the electric field where E = E/4, as can be seen by taking,
for example, the field in the z direction:

Ez (r , z ) =

2Qi z E(k)
1 2Qi z E(k)
1
=
= Ez (r, z),
3
2
3
2
π R 1−k
4 π R 1−k
4

(2.33)

where we have performed the substitutions z → 2z, R → 2R. A similar calculation can be
done for the field in the radial direction. This result can be explained intuitively by noting
that the units of the electric field are [charge]/[length]2 . If the length is doubled, the electric
field scales as one fourth of the original value.
In practice for the simulations presented here, this means that if the radius of the surface
charge increases as r0 = αr0 , it is sufficient to increase the coordinates corresponding to the
mesh points as r = αr, z = αz and scale the electric field at each point as Er,z = Er,z /α2 ,
where Er,z is a value calculated and tabulated before each simulation run. This method has
the effect of increasing the mesh spacing as the pulse expands, with the advantage that the
values for the positive field at a point need to be calculated only once, saving computational
time and resources.

2.5

Pulse time evolution

We now turn our attention to the results of these N-particle simulations, starting from an
analysis of the time dependent evolution of a single pulse followed by a comparison of the
properties of the emitted pulses for varying extraction conditions.
We begin by presenting, in Fig. 2.3, the temporal dependence of various key pulse

22

parameters for Ne0 = 7 · 106 generated electrons with extraction fields of Fa =0.32 MV/m and
Fa =10 MV/m. In panel (a) of the figure, the number of electrons Ne is plotted as a function
of time on a log-log scale, showing the increase in Ne during the initial photoemission phase
lasting up to 0.1 ps. During the time evolution of the pulse, the electrons can recombine
with the photocathode surface, leading to a decrease in Ne at later times. This effect plays
a significant role for low extraction fields since the pulse remains closer to the surface, and
is absent in the high field case (dashed line in the figure). Figure 2.3(b) shows the pulse
center of mass distance from the surface for the two extraction fields considered compared to
a free particle moving in the corresponding electric field [z = 1/2(eFa /m)t2 ]. The center of
mass evolution for Fa =10 MV/m is free particle like, while at lower values of Fa , the positive
surface field hinders the pulse propagation.
Next, we look at the pulse size and normalized emittance [panels (c) and (d)], two parameters that are relevant for our imaging and diffraction applications. Looking at the transverse
components of these quantities [black lines in panels (c), (d)], we note that they remain close
to the initial values for both data sets presented, with a weak dependence on extraction field,
seen mostly in the increase of the pulse width σx for t >50 ps and Fa =10 MV/m due to the
weaker effect of the positive surface charge, which for Fa =0.32 MV/m keeps the electrons
focused towards the central region of the pulse. On the other hand, the longitudinal components of both pulse size and emittance [red lines in panels (c), (d)] show a very rapid increase
during the first 20 ps of the simulation, due to the strong Coulomb repulsion between the
electrons extracted from the photocathode surface and even after this initial expansion stage,
both longitudinal components keep increasing for the duration of the simulation.
To better understand these features we look at the time evolution of the pulse in the
x − z plane, shown in Fig. 2.4 with a color map of the charge distribution in the rest
23

(b) 105

(a)
10

104

106

103

zcm (µm)

Ne

7

105

2

10

101
4

10

103
0.001 0.01

0.1
1
t(ps)

10

Fa = 0.32MV/m
Fa = 10MV/m
free particle

100

Fa = 0.32MV/m
Fa = 10MV/m

10-1

100

(c)

0

20

40

60
t(ps)

80

100 120

(d)
10-1

101

Fa = 0.32MV/m
σx
σz

100

Fa = 10MV/m
σx
σz

-1

10

10-2

0

20

40

60
t(ps)

80

10-3

Fa = 0.32MV/m
εx
εz

10-5

Fa = 10MV/m
εx
εz

εi (µm)

σi (µm)

102

100 120

10-7

0

20

40

60
t(ps)

80

100 120

Figure 2.3: Time evolution of key pulse parameters, for Ne0 = 7 · 106 , extraction fields of
Fa =0.32 MV/m (continuous lines) and Fa =10 MV/m (dashed lines): (a) number of electrons
Ne (note the logarithmic scale of both axes). (b) Pulse center of mass distance from the
surface and corresponding free particle time evolution [z = 1/2(eFa /m)t2 , dotted lines]. (c)
Transverse (black) and longitudinal (red) pulse size. (d) Transverse (black) and longitudinal
(red) normalized emittance.

24

frame of the bunch overlaid with arrows representing the average electron velocity. For short
timescales we observe that the spatial profile of the pulse retains the Gaussian distribution of
the generating laser pulse [Fig. 2.4(a)] and the velocities show a turbulent flow related to the
initial thermal distribution from which they were generated and independent of extraction
field, since the repulsion due to the space charge has not had time to play a big role. At
later times the pulse begins to move away from the surface, with a significant residual charge
left behind for low extraction fields [Fig. 2.4(b)]. The final profile of the pulse as t=120 ps,
depends strongly on the extraction field. Comparing panels 2.4(c) and 2.4(d), corresponding
respectively to Fa =0.32 MV/m and Fa =10 MV/m, we see that in the latter picture the flow
has become fully laminar and the pulse has evolved into the typical pancake shape, while
for the lower extraction field the pulse retains some turbulent flow and does not fully detach
from the surface.
Figure 2.4(e) presents the radial pulse profile at t=120 ps for three different extraction fields fitted to a combination of Gaussian [G = a exp(r2 /2σ 2 )] and Ellipsoidal [E =
a 1 − (r/R)2 ] functions. The laminar flow observed in the charge distribution plots [panel
(d)] can be seen to correspond to the formation of an Ellipsoidal profile, which is advantageous due to its linear internal force field that leads to a non increasing emittance [53, 64].
As a consequence even though the pulse will expand and change its dimensionality ratio, this
expansion can be fully corrected with appropriate beam optics such as compression [65] and
aberration corrections [40]. On the other hand at Fa =0.32 and 1 MV/m, the flow is not fully
laminar; this component can be described by a Gaussian with extended tails representative
of the thermal-like population (accounting for, respectively 35% and 27% of the total number
of electrons), which coexists with the elliptical component.
Overall the emergence of a 2D ellipsoidal profile from the complex interplay of image- and
25

Figure 2.4: (a)-(d) Color map of the charge distribution in the rest frame of the bunch
projected onto the x-z plane, overlaid with arrows representing the average electron velocity.
Panel (a) shows the initial stage of photoemission, which is independent of extraction field.
Panel (b), (c) show the intermediate and final stages for weak (Fa =0.32 MV/m) extraction
field at times t=60 ps and t=120 ps. Panel (d) corresponds to the final stage for the strong
field case (Fa =10 MV/m). Note that the number of emitted electrons is the same in all
cases, Ne0 = 7 · 106 , while the final number of emitted electrons depends on the extraction
field. (e) Transverse charge density profiles under different extraction fields (Fa =0.32, 1 and
10 MV/m) fitted with a combination of Gaussian (G) and Ellipsoidal (E) functions [26].
26

space charge effects, is somewhat surprising since formation of an ellipsoidal beam driven by
self fields is expected only at zero surface charge limit, i.e. σ0 /

0

Fa , where

0

is vacuum

permittivity and σ0 = 3Ne e/(2πr2 ) is the surface charge density. Inserting the parameters
used for the simulations in Fig. 2.4 into the equation gives σ0 =5.35 × 10−5 C/m2 and σ0 /

0

=

6.04 MV/m, showing that the extraction fields used (Fa ) are comparable to this value. With
the highest number of particles considered, Ne0 = 2 · 108 , σ0 /

0

= 172.55 MV/m, and this

condition is definitely not satisfied. This is in contrast with the situation described elsewhere
in the literature such as in the work by Luiten et al. [53, 64] where the extraction field Fa =
100 MV/m and σ0 /

0

= 5 MV/m, in which case the bunch evolution is solely dependent on

its self fields as image- and space-charge effects are negligible.
We conclude this section with a comparison of the simulated data to the values measured
in the shadow imaging experiments [66] (Fig. 2.5). Panel (a) shows the raw experimental
images, while panels (b)-(d) compare the longitudinal pulse profile in experiment (red circles) and simulation (black line) at various times. The shaded box indicates the region in
the proximity of the surface, which cannot be accurately imaged with the shadow imaging
technique. The experimental data is shifted horizontally to fix the position of the surface
and scaled vertically to match the units used in the simulation. This fitting is done for the
data set at t=50 ps (not shown) and fixed for the remainder of the plots. Given the few
parameters used in the simulation, we see a remarkable agreement with the measured pulse
profile over a time period of t

50 ps. This provides a validation that the simulation method

used correctly describes the experimental system.

27

Figure 2.5: Comparison to experimental data. Panel (a) shows the raw experimental images,
while panels (b)-(d) show a comparison of the longitudinal pulse profile at various times in
experiment (red circles) and simulation (black line), corresponding to an extraction field of
Fa = 0.32 MV/m and Ne0 = 108 [26].

2.6

Optimal pulse generation

In this section we investigate the optimal conditions for pulse generation in terms of laser
fluence (or, equivalently, initial number of electrons, Ne0 ) and extraction field Fa for different
experimental realizations. We start the discussion by presenting our results for photoemission
using a Gaussian laser pulse from a flat surface (modeled with a radially expanding positive
surface charge), followed by a comparison with the pinned case (stationary positive charge
on the surface). To understand the effect of the spatial profile of the laser, we also discuss
our results for the radially Uniform (top-hat) and Elliptical profiles.

28

2.6.1

Gaussian pulse

In order to use the generated electron bunch for single-shot imaging, a high enough number of
electrons is necessary, so one of the first quantities of interest is the final number of electrons
at time t=120 ps. This is shown in Fig. 2.6 as a function of the initial number of electrons, Ne0
for extraction fields Fa = 0.1 MV/m to 40 MV/m. For a low number of generated electrons
(Ne < 106 ), all of the particles escape the surface, creating the linear relationship seen in the
first part of the figure. As the initial number of electrons is increased, the space charge effect
also increases due to both the electrons and the image charge on the surface. This creates the
so-called virtual cathode (VC) limit in which the negative charge of the electrons emitted at
earlier times, combined with the attractive surface field, hinders further emission of particles
and effectively limits the final number of electrons extracted, a feature that can be seen in
the flattening of the curves shown. Increasing the extraction field shifts the onset of this
phenomenon enabling the emission of a higher number of particles from the photocathode.
The dependence of the onset of the virtual cathode effect on extraction field is investigated
by calculating the critical number of electrons from a linear fit of the two regimes described.
The resulting data is plotted in the inset of Fig. 2.6 as a function of extraction field Fa . If
we approximate the electron pulse as a sheet of charge and define the onset of the virtual
cathode regime to be when the electric field generated by the sheet of emitted electrons is
equal to the extraction field Fa [67], we can express the critical current required as:

Jcrit =

0 Fa /τp ,

(2.34)

where τp is the pulse duration. Having assumed a thin sheet of charge, we can write σ = Jτp ,

29

where σ = Necrit /(πr2 ) is the electron charge density in the xy plane so that

σ=

0 Fa .

(2.35)

The fit of our data to this equation is also shown in the inset of Fig. 2.6, with a numerical
prefactor of α = 1.65 so that σ = α 0 Fa (α = 1 in the ideal case).

8

10

Ne

crit

108
107

106

107

0.1

1

10

Neemit

Fa (MV/m)

106

∝ (Ne0)1
∝ (Ne0)1/3

Fa (MV/m)
0.1
0.2
0.32
1
10
40

105
105

106

107
Ne

108

0

Figure 2.6: Number of electrons emitted at t=120 ps as a function of number of generated
electrons for various extraction fields. The lines serve as a guide to the eye for the scaling
Neemit ∝ Ne0 and Neemit ∝ (Ne0 )1/3 . The data shows evidence of virtual cathode formation.
Inset: Critical number of electrons for virtual cathode formation as a function of extraction
field Fa . The dashed line shows the critical number of electrons as a function of Fa obtained
by using Eq. 2.35 [26].

The onset of the VC limit also acts on the longitudinal pulse width, σz (Fig. 2.7). For
pulses in the sub-VC regime, σz increases linearly with the number of particles, while in
the VC regime the width of the pulse does not show a strong dependence on the number

30

of electrons emitted. Our simulations show that the trend of σz ∝ (Neemit )1/2 seen in the
experimental data (right panel of Fig. 2.7), is a result of the linear growth of σz in the
sub-VC regime combined with its reduced increase as the VC regime sets in.

Figure 2.7: Longitudinal pulse width, σz , as a function of number of electrons at time
t=120 ps, for various extraction fields Fa from the simulation (left panel) and compared to
experimental data (right panel). The lines serve as a guide to the eye for the cases σz ∝ Neemit
and σz ∝ (Neemit )1/2 [26].

2.6.2

Effect of image charge pinning

We now focus on the effect that pinning of the image charge on the surface has on the onset
of the VC limit and on the final pulse parameters, as this will provide information on the
optimal photocathode geometry one should employ. By looking at the charge distribution
and transverse profiles of individual pulses (Fig. 2.8), we observe that in the pinned case
the positive image charge field provides a focusing effect, increasing the fraction of electrons
at the center of the pulse and consequently reducing the fraction of electrons in the lateral
region. The particles in the central region are also pulled back towards the surface by the
strong electric field, creating a channel of electrons that prevents the pulse from becoming
31

fully pancake-like even at t=120 ps.
Fa=10 MV/m

Fa=10 MV/m - Pinned

(a)

(b)

200

200
x (µm)

400

x (µm)

400

0

0

-200

-200

-400

-400

-600

-300

0
z (µm)

Fa = 0.32 MV/m

300

-600

Fa = 1.0 MV/m

-300

0
z (µm)

300

Fa = 10 MV/m

(c)
Electron fraction (a.u.)

Expanding
Pinned

-400 -200

0 200 400 -400 -200 0 200 400 -400 -200 0 200 400
x (µm)
x (µm)
x (µm)

Figure 2.8: (a),(b) Color map of the charge distribution in the rest frame of the bunch
projected onto the x-z plane, overlaid with arrows representing the average electron velocity.
Panel (a) shows the distribution at t=120 ps for extraction field Fa = 10 MV/m and number
of emitted electrons Ne0 = 7 · 106 . Panel (b) shows the charge distribution under the same
conditions but with the image charge pinned to a small area on the photocathode surface. (c)
Transverse charge density profiles under different extraction fields (Fa =0.32, 1 and 10 MV/m)
showing the central focusing peak caused by the pinned image charge.

Image charge pinning has very little effect on the longitudinal pulse width σz or on the
number of electrons emitted from the surface before the onset of the VC (Fig. 2.9). Upon
reaching the VC limit, though, the number of electrons in the pinned case is reduced by
≈ 5% compared to the non-pinned case, while maintaining the same linear dependence on

32

Ne0 seen in Fig. 2.6.
(a)

(b)

107

Pinned
Fa (MV/m)
0.1
0.2
0.32
1
10
40

σz (µm)

Neemit

108

Pinned
Fa (MV/m)
0.1
0.2
0.32
1

102

106
Expanding

Expanding
105

1

105

106

107
Ne

10

108

105

0

106
Ne emitted

107

Figure 2.9: (a) Number of electrons emitted at t=120 ps as a function of number of generated
electrons for various extraction fields in the presence of a pinned image charge. The red lines
show the data corresponding to the expanding image charge. (b) Longitudinal pulse width,
σz , as a function of number of electrons at time t=120 ps, for various extraction fields Fa with
a pinned image charge. The red lines, again, show the data corresponding to the expanding
image charge.

2.6.3

Effect of transverse laser profile

We now turn our attention to the effect of the laser pulse shape on the properties of the
photoemitted electron pulse and start by observing the temporal evolution of an electron
bunch generated from a radially uniform laser pulse (Fig. 2.10), where the laser intensity
I(r) ∝ const. Also shown for comparison is the time evolution of a similarly generated pulse
in the absence of the surface image charge. Similarly to the pinned case, the positive field on
the surface attracts the electrons and causes a focusing of the electron density towards the
center of the pulse where the field is strongest [Fig. 2.10(e)]. This effect is seen to be more
pronounced for low extraction fields Fa , as the pulse remains close to the surface throughout
the simulation.
In the case of an elliptical laser pulse we have that the intensity I(r) ∝
33

1 − (r/R)2 ,

(a)

(b) 600
200

t=0.63 ps

t=120 ps
300
x (µm)

x (µm)

100
0
-100

0
-300

-200
-0.05

0
z (µm)

0.05

-600
-900 -600 -300
0
z (µm)

0.1

(c)

(d)
200

t=0.63 ps

600

t=120 ps

300
x (µm)

100
x (µm)

600

300

0

0

-100

-300

-200

-600
-0.1

-0.05

0
z (µm)

Fa = 0.32 MV/m

0.05

0.1

-600

Fa = 1.0 MV/m

-300

0
z (µm)

300

600

Fa = 10.0 MV/m

Electron fraction (a.u.)

(e)
With IC
No IC

-600 -300

With IC
No IC

With IC
No IC

0 300 600 -600 -300 0 300 600 -600 -300 0 300 600
x (µm)
x (µm)
x (µm)

Figure 2.10: (a)-(d) Color map of the charge distribution in the rest frame of the bunch
projected onto the x-z plane, overlaid with arrows representing the average electron velocity
for the top-hat (uniform) pulse. Panels (a) and (b) show the initial and final stages of
photoemission in the presence of an image charge on the surface. Panels (c), (d) correspond
to the case without image charge. Note that the number of emitted electrons is the same in
all cases, Ne0 = 1 · 107 , with an extraction field Fa = 1 MV/m . (e) Transverse charge density
profiles under different extraction fields (Fa =0.32, 1 and 10 MV/m), showing the effect of
the image charge.

34

where R corresponds to the spatial width of the laser. Figure 2.11 shows the phase space
plots and radial distributions corresponding to this case both with and without the image
charge. Again we observe a dramatic change in the overall pulse shape that is a consequence
of the drastic disruption of the linear self fields otherwise present in the elliptical pulse [53].
Projecting the charge distribution onto the x-axis [Fig.2.11(e)], show that this effect results
in the formation of a pronounced central peak. A fit of the data to Gaussian and elliptical
distributions shows an increasing weight of the Gaussian pulse component at low extraction
fields (67% at Fa = 0.32M V /m and 37% at Fa = 10M V /m) for which the perturbation to
the pulse profile is more pronounced as the electrons remain closer to the surface throughout
the simulation.

2.7

Guidelines for generation of high brightness electron beams

Next we wish to compare the pulse properties obtained under these different conditions of
pulse generation. As our goal is to be able to do time resolved imaging, both the transverse
and longitudinal emittances are of importance since they determine, respectively, the spatial
and temporal resolutions. We also compare the emittance, coherence length and energy
spread obtained in each case.
Figure 2.12 shows the emittances corresponding to the Gaussian laser pulse both with
expanding and pinned surface charges. Our data shows a strong increase of the transverse
emittance

x

[panel (a)] upon reaching Necrit , due to the onset of the VC regime. From the

inset, which shows the temporal evolution of
state value of

x

x

at Fa =0.32, 1 MV/m, we see that a steady

is reached after approximately 40 ps and is correlated with the region where
35

(a)

(b) 600
200

t=0.63 ps

t=120 ps
300
x (µm)

x (µm)

100
0
-100

0
-300

-200
-0.05

0
z (µm)

0.05

-600
-900 -600 -300
0
z (µm)

0.1

(c)

(d)
200

No IC - t=0.63 ps

600

No IC - t=120 ps

300
x (µm)

100
x (µm)

600

300

0

0

-100

-300

-200

-600
-0.1

-0.05

0
z (µm)

Fa = 0.32 MV/m

0.05

0.1

-600

Fa = 1.0 MV/m

-300

0
z (µm)

300

600

Fa = 10.0 MV/m

Electron fraction (a.u.)

(e)
With IC
No IC

-600 -300

With IC
No IC

With IC
No IC

0 300 600 -600 -300 0 300 600 -600 -300 0 300 600
x (µm)
x (µm)
x (µm)

Figure 2.11: (a)-(d) Color map of the charge distribution in the rest frame of the bunch
projected onto the x-z plane, overlaid with arrows representing the average electron velocity
for the elliptically shaped pulse. Panels (a) and (b) show the initial and final stages of
photoemission in the presence of an image charge on the surface. Panels (c), (d) correspond
to the case without image charge. Note that the number of emitted electrons is the same in
all cases, Ne0 = 1 · 107 , with an extraction field Fa = 1 MV/m. (e) Transverse charge density
profiles under different extraction fields (Fa =0.32, 1 and 10 MV/m), showing the effect of
the image charge.

36

the previously discussed 2D elliptical pulse structure dominates. In contrast,

z

[Fig. 2.12(b)]

shows no change in behavior as the VC regime is reached and is also rather insensitive to
the extraction field. This is an indicator that in our pancake-like regime, the transverse
emittance is driven by the strong nonlinear fields while the longitudinal component has a
universal dependence on Neemit . This is particularly evident in the inset, showing

z (t)

for

three very different external field values and Neemit = 107 . Pinning of the surface charge has
the effect of increasing the transverse emittance significantly once the virtual cathode limit
is reached; the longitudinal emittance is also similarly affected.
Changing the shape of the laser profile used to emit the electrons also has a significant
impact on the normalized emittance (Fig. 2.13). For Fa = 1 MV/m (black symbols and lines),
the elliptical shape (circles) is favorable below the VC limit as it has a smaller transverse
emittance [panel (a)] compared to the Gaussian and uniform (top-hat) profiles. Previous
literature [53, 64], showed this to be the case for negligible image charge fields due to the
internal linear self-forces that limit the increase of the emittance. It is non-trivial that this
would still be the case even in the presence of the positive field on the surface, which disrupts
the linearity of the electron-electron forces and drastically changes the resulting pulse profile,
as we saw in Fig. 2.11. The results in Fig. 2.13 show that once the VC limit is reached,
at about Neemit

1 · 107 for Fa = 1 MV/m, both the elliptical and uniform profiles present

a sharp scaling with the number of electrons, while the Gaussian profile retains a more
favorable scaling, with a slope for the emittance increase that is 1/3 that of the elliptical
case. The origin of this scaling can be understood by looking at the inset in Fig. 2.13, which
presents the number of emitted electrons Neemit as a function of Ne0 , the initial number of
electrons. In the elliptical and uniform case, the VC effect completely blocks the increase
in Neemit , such that for Ne0 > 107 , the curves shown in the inset are flat. With a Gaussian
37

(a)

1
εx (µm)

100
10-1

1
0.32
0.32 P.

-2

10

-3

10

20

40 60 80 100 120
time (ps)

εx (µm)

0

Fa (MV/m)
0.1
0.32
1.0
10
0.32 - Pinned

0.1

(b)

1

εz (µm)

100
10-2

1
0.32
0.32 P.

10-4
10-6

εz (µm)

0

0.1

20

40 60 80 100 120
time (ps)

Fa (MV/m)
0.1
0.32
1.0
10
0.32 - Pinned

0.01

0.001
105

106

107
Ne

108

emit

Figure 2.12: (a) Transverse normalized emittance x and (b) longitudinal normalized emittance z as a function of number of emitted electrons Neemit and extraction field Fa at 120 ps
for a Gaussian pulse. Insets: time dependence of x and z for three selected cases at
Neemit = 107 [26].
laser profile the emission of electrons is, instead, only hindered allowing an extraction of a
higher number of particles under the same conditions. An increase in the extraction field
Fa has the effect of shifting the onset of the VC regime to higher Neemit , thus increasing the
38

regime before the VC is reached in which the elliptical pulse is optimal. The longitudinal
emittance

z

[Fig. 2.13(b)] scales linearly with the number of electrons for all laser profiles

considered, with a weak dependence on the extraction field. The signature of the onset of
the VC regime is also seen here in the sudden increase of
for Neemit

z

by about one order of magnitude

1 · 107 .

The influence of the image charge on the emittance can be understood by comparing the
time evolution of the 6D normalized emittance (

x

·

y

· z ) as a function of laser pulse profile

and presence of image charge, shown in Fig. 2.14 both below and above the VC limit. Below
the virtual cathode limit [Fig. 2.14(a)], the image charge acts as a weak perturbation on the
pulse dynamics with some effect only during the initial stages (t <60 ps) of pulse formation.
In contrast for extraction conditions above the virtual cathode limit [panel (b)], the strong
positive field on the surface has a small effect on the emittance of the Gaussian pulse profile
(black continuous and dashed lines) but is a significant perturbation in the case of both the
elliptical (red) and uniform (green) cases.
From this we can conclude that the optimal extraction conditions depend on a delicate
balance between the image charge, which disrupts the favorable internal self-fields of both the
elliptical and uniform cases, the number of electrons emitted and the extraction field, which
controls the onset of the virtual cathode limit. Below this limit the image charge has a weak
effect with the pulse properties primarily controlled by the laser parameters (laser profile
and fluence), while above the VC limit the image charge strongly perturbs the expansion of
the pulse with a nonuniform dependence on the pulse profile.
A further check of these considerations is shown in Fig. 2.15 where the 6D emittance is
plotted as a function of time for the Gaussian and elliptical cases using different models for
the image charge: expanding radially, pinned, and placed behind the surface at a distance
39

1

107

εx (µm)

Ne

emit

(a)

0.1

Gauss.
Ellipt.
Unif.

106

106

107

Ne

108
0

Fa = 1.0 MV/m
Gaussian
Elliptical
Uniform
Fa = 10 MV/m
Gaussian
Elliptical
Uniform

(b) 0.01

0.1

Fa = 1.0 MV/m
Gaussian
Elliptical
Uniform

εz (µm)

1

Fa = 10 MV/m
Gaussian
Elliptical
Uniform

0.01

0.001

105

106

107
Ne

108

emit

Figure 2.13: (a) Transverse normalized rms emittance x and (b) longitudinal normalized
rms emittance z dependence on number of emitted electrons Neemit and extraction field Fa
for varying shapes of the exciting laser pulse at t = 120ps. Inset: number of emitted electrons
Neemit at t=120 ps as a function of Ne0 , the initial number of electrons for Fa = 1 MV/m and
varying laser profile.
z = −zCM , where zCM is the distance of the center of mass of the electron pulse. Below the
virtual cathode limit, the emittance values are robust with respect to the different choices

40

(a)

-6

10

εx εy εz (µm3)

10-7

-8

10

Below VC
G
G0
E
E0
U
U0

10-9
(b)

10-2

εx εy εz (µm3)

10-3
10-4
10-5
Above VC
G
G0
E
E0
U
U0

10-6
10-7
10-8

0

20

40

60

80

100

120

time (ps)

Figure 2.14: Time dependence of the 6D normalized emittance for Fa = 1 MV/m both (a)
below and (b) above the virtual cathode limit. The data shown corresponds to the Gaussian
(G), Elliptical (E) and Uniform (U) cases, with the subscript 0 indicating the absence of the
image charge (G0, E0, U0 respectively).
in treating the image charge. On the other hand in the nonlinear regime shown by the pulse
above the VC limit, the choice of the model for the image charge can have the effect of
varying the final value of emittance by about 2 orders of magnitude for the elliptical case,
41

and by about one for the Gaussian profile. Conventionally UEM systems always operate
below the VC regime where the image charge field does not have a critical effect, but this
raises an important issue that needs to be addressed in pushing to higher brightness electron
pulses.
To conclude our discussion and extract useful guidelines for electron beam generation, we
compare the brightness (Fig. 2.16), coherence length and energy spread (Fig. 2.17) for pulses
generated under the different conditions described in the previous sections. Since pinning
of the image charge on the surface has a negligible effect below the VC and an unfavorable
scaling after the VC, it has been omitted for clarity from the discussion that follows.
The first parameters we wish to look at are the 4D and 6D brightness defined as B4D =
Ne /(

x

·

y)

and B6D = Ne /(

x

·

y

·

z ),

reported in Fig. 2.16. If temporal resolution is not

required, the appropriate figure of merit is the 4D brightness [panel (a)], where we observe
three distinct regimes for a given extraction field Fa . At low pulse currents (Ne < 104 for
Fa = 1 MV/m), it is independent on the transverse pulse shape. For intermediate values
of Ne (104 < Ne < 3 · 107 for Fa = 1 MV/m), an elliptical spatial profile offers the best
solution. At high numbers of electrons (Ne > 3 · 107 for Fa = 1 MV/m), the Gaussian
pulse offers a better scaling than the Elliptical or Uniform but is associated with a very low
peak brightness. Increasing the extraction field solves this issue by shifting the boundaries
between the regimes discussed to higher values of Ne . On the contrary, if temporal resolution
is a concern, it is necessary to include the longitudinal pulse properties in the discussion and
therefore the 6D brightness is the relevant figure. From our results in Fig. 2.16(b), we
conclude that in this case an elliptical transverse profile offers the optimal choice for the
whole parameter range considered.
Other key quantities for electron source design are the coherence length (Eq. 2.1) and
42

(a)

-6

10

-7

εx εy εz (µm3)

10

Below VC
Gaussian
G. Pinned
G. Behind Surf.
Elliptical
E. Pinned
E. Behind Surf.

-8

10

(b)

10-9
10-2

εx εy εz (µm3)

10-3
10-4
10-5
Above VC
Gaussian
G. Pinned
G. Behind Surf.
Elliptical
E. Pinned
E. Behind Surf.

-6

10

10-7
10-8

0

20

40

60

80

100

120

time (ps)

Figure 2.15: Time dependence of the 6D normalized emittance for Fa = 1 MV/m both
(a) below and (b) above the virtual cathode limit. The data shown corresponds to the
Gaussian (G) and Elliptical (E) cases, for different ways of considering the image charge
field: expanding, pinned to the surface or mirroring the electron pulse distance from the
surface.
the energy spread (Eq. 2.3), These two quantities are plotted in Fig. 2.17 as a function of
emitted electron number and photoemitting laser pulse shape. The Elliptical profile presents

43

4D Brightness (x108 µm-2)

(a) 100
Pancake:
Gaussian
Fa (MV/m)
1.0
10
Elliptical
Fa (MV/m)
1.0
10

10

1

Uniform:
Fa (MV/m)
1.0
10

0.1

0.01
(b)

6D Brightness (x109 µm-3)

1000

Pancake:
Gaussian
Fa (MV/m)
1.0
10
Elliptical
Fa (MV/m)
1.0
10

100

10

Uniform:
Fa (MV/m)
1.0
10

1

0.1

0.01

103

104

105

106
Ne

107

108

emit

Figure 2.16: (a) 4D brightness and (b) 6D brightness as a function of number of electrons
emitted, Neemit for the Gaussian, Uniform and Elliptical cases. The data shown corresponds
to Fa = 1 MV/m (closed circles) and Fa = 10 MV/m (open circles).
the highest coherence length compared to the other two data sets even after reaching the
VC limit. Increasing the extraction field offers the advantage of pushing the onset of the VC
limit to higher Ne , increasing the range over which the coherence length is constant. In terms

44

of the energy spread, the Elliptical and Uniform profiles are both equally favorable, having
a lower ∆E than the Gaussian case and, overall, this parameter is relatively insensitive to
the extraction field.

10

Coherence length (nm)

(a)

(b)

Pancake:
Gaussian
Fa (MV/m)
1.0
10
Elliptical
Fa (MV/m)
1.0
10
Uniform:
Fa (MV/m)
1.0
10

1

105

Pancake:
Gaussian
Fa (MV/m)
1.0
10
Elliptical
Fa (MV/m)
1.0
10

104

∆E (eV)

103
102

Uniform:
Fa (MV/m)
1.0
10

101
100
10-1

103

104

105

106
Ne

107

108

emit

Figure 2.17: Coherence length (a) and energy spread (b) as a function of number of electrons
emitted, Neemit for the Gaussian, Uniform and Elliptical profiles. The data shown corresponds
to Fa = 1 MV/m (closed circles) and Fa = 10 MV/m (open circles).

45

Chapter 3
Analytic Gaussian Model
The photoemission process previously described is only the first step in the operation of an
Ultrafast Electron Microscope (UEM) system. We now use the knowledge acquired from the
simulations of pulse formation as an input for calculations involving the whole microscope
column, which was introduced in Chapter 1. To do so, we use the so-called Analytic Gaussian
model introduced by Michalik and Sipe [37,38] and refined by Berger and Schroeder [39]. We
begin the discussion with an overview of the derivation of the equations for this mean-field
model and the proceed to show how we added the contributions due to the photoemission
gun and the relativistic motion of the electrons to this formalism to make it applicable for
our system. We then explain the procedure used for optimizing the lens and RF cavity
strengths and analyze the results obtained, taking into account the limits inherent in this
model.

46

3.1

Derivation of the Analytic Gaussian equations

This section elucidates the derivation of the mean field Analytic Gaussian differential equations, following the work in [37]. We start with the propagation of an N-electron pulse in
free space, and then add in the effects of lenses, radio frequency (RF) cavities, photoemission
gun and relativistic motion.

3.1.1

N-electron pulse in the absence of external forces

The equations of motion for a classical system of N electrons are given by

p (t)
dri (t)
= i ,
dt
m
∂V [ri (t) − rj (t)]
dpi (t)
=−
,
dt
∂ri (t)
j=i

(3.1)
(3.2)

where ri are the positions of particles i = 1, ..., N , pi are the momenta, m is the electron
mass and V (r) is the Coulomb potential

V (r) =

with e the electron charge and

0

e2
,
4π 0 |r|

(3.3)

the vacuum permittivity. We now introduce a discrete

distribution to describe the particles in the system, which we formally write as
N

δ [ri (t) − r] δ [pi (t) − p] ,

fD (r, p; t) =
i=1

47

(3.4)

where fD is normalized so that

drdpfD (r, p; t) = N.

(3.5)

The Boltzmann equation for this distribution can then be written as

∂fD (r, p; t)
∂fD (r, p; t)
=−
∂t
∂r
∂fD (r, p; t)
=−
∂r

∂r ∂fD (r, p; t) ∂p
−
·
,
∂t
∂p
∂t
p
∂fD (r, p; t)
·
−
· F,
m
∂p
·

(3.6)
(3.7)

where Eqs. 3.1 and 3.2 have been used. We rewrite the force term

F =−
j=i

∂
∂r
∂
=−
∂r

=−

∂V [ri (t) − rj (t)]
,
∂ri (t)

(3.8)

dr ρ(r )V (r − r ),

(3.9)

dr dp fD (r , p ; t)V (r − r ),

(3.10)

where ρ(r ) is the local density equal to

ρ(r ) =

dp fD (r , p ; t),

(3.11)

and then substitute Eq. 3.10 into Eq. 3.7 to obtain

∂fD (r, p; t)
∂fD (r, p; t) p
∂fD (r, p; t) ∂
=−
·
+
·
∂t
∂r
m
∂p
∂r

dr dp fD (r , p ; t)V (r − r ). (3.12)

The discrete distribution fD depends on the specific initial positions and momenta of

48

the particles, so next we take an ensemble average over a set of initial conditions, which we
indicate as f (r, p; t) = fD (r, p; t) . Performing this average, Eq. 3.12 becomes

∂f (r, p; t)
∂f (r, p; t) p
=−
·
+
∂t
∂r
m

∂fD (r, p; t) ∂
·
∂p
∂r

dr dp fD (r , p ; t)V (r − r ) . (3.13)

As is typical in this kind of problem, it is impossible to derive a closed form expression
without solving for higher order moments such as fD (r, p; t)fD (r , p ; t) . By using a mean
field approach, we approximate fD (r, p; t)fD (r , p ; t)

fD (r, p; t) fD (r , p ; t) , so that

Eq. 3.13 becomes

∂f (r, p; t)
∂f (r, p; t) p
∂f (r, p; t) ∂
=−
·
+
·
∂t
∂r
m
∂p
∂r

dr dp f (r , p ; t)V (r − r ).

(3.14)

Even with this approximation, it is not possible to solve Eq. 3.14 analytically, so we
choose to parameterize the distribution f (r, p; t) with some time-dependent parameters for
which we can derive ordinary differential equations. Using the cylindrical symmetry present
in the UEM system and assuming for simplicity a Gaussian distribution for both the transverse and longitudinal components, we write

f (r, p; t) = C exp −

x2 + y 2
z2
[px − (γT /σT )x]2 + [py − (γT /σT )y]2 [pz − (γz /σz )z]2
−
−
−
,
2σT
2σz
2ηT
2ηz
(3.15)

where σT,z , ηT,z and γT,z are time dependent parameters that describe, respectively, the
spatial width, kinetic energy spread and chirp of the electron pulse (see Fig. 3.1). Note
that the units for the parameters are σT,z units of length squared, ηT,z units of momentum
squared and γT,z units of action (length x momentum). C is a time dependent normalization

49

constant
C(t) =

N
(2π)3

1
2 2
σT ηT σz ηz

1/2

.

(3.16)

γi
σi

pi

xi

ηi1/2

σi1/2
Figure 3.1: Schematic representation of the pulse parameters used in the Analytic Gaussian
model in phase space.
We can rewrite Eq. 3.15 in a more general form as a multivariate normal distribution

f (u; t) =

N
(2π)k |χ|

1
exp − uT χ−1 u ,
2

(3.17)

where u = (x, px , y, py , z, pz ), k is the dimension of u (k = 6 in our case), χ is the covariance matrix and |χ| is its determinant. The distribution is multiplied by N to obey the
normalization constraint in Eq. 3.5. Comparing Eqs. 3.15 and 3.17, we see that

χ−1



 aT 0 0


=  0 aT 0


0 0 az
50




,



(3.18)

with


ai = 

1
σi

+

γi2
ηi σi2

γi
ηi σi

γi
ηi σi
1
ηi



,

(3.19)

where i = T or z. Calculating the covariance matrix χ is trivially done by inverting the
matrix in Eq. 3.18 giving




a−1
T

0




χ= 0


0

a−1
T
0

0 


,
0 


a−1
z

(3.20)

where a−1
represents the inverse of matrix ai :
i




γi

 σi
ai−1 = 
,
2
γi
γi ηi + σi

(3.21)

where, again, i = T or z.
Using the definition of the covariance matrix [68], we know that

χi,j =< ui uj >=

1
N

du ui uj f (u; t),

(3.22)

where we use < ... > to represent an expectation value over the ui , uj components of the
vector u = (x, px , y, py , z, pz ). We now take a time derivative of Eq. 3.22:

∂
1
χi,j =
∂t
N

du ui uj

∂f (u; t)
,
∂t

(3.23)

and substitute ∂f (u; t)/∂t with the mean-field equation 3.14 to obtain differential equations

51

for each matrix component χi,j

∂
χi,j = Kijf low + Kijf orce ,
∂t

(3.24)

where

Kijf low = −
Kijf orce =

1
N

1
N

p ∂f (u; t)
·
,
m
∂r
∂f (u; t) ∂Φ
·
,
du ui uj
∂p
∂r
du ui uj

(3.25)
(3.26)

with the notation u = (x, px , y, py , z, pz ), r = (x, y, z), p = (px , py , pz ) and

Φ=

du f (u ; t)V (r − r ).

(3.27)

Note that due to the symmetry of the problem, the terms corresponding to the x and y
coordinates are equal, meaning that χ11 = χ33 , χ22 = χ44 and χ12 = χ21 = χ34 = χ43 .
Similar relations hold for the Kij terms. For this reason in the following discussion we will
only consider components in the x direction.
The K f low terms describe the geometric motion of the electrons and are easily evaluated
by substituting the explicit form of f (u; t) (Eq. 3.15) and solving the resulting Gaussian
f low
f low
f low
integrals. It is easy to see that the K22
, K44
, K66
elements vanish and that the

52

remaining terms are

2
γT ,
m
1
γ2
=
ηT + T
m
σT

f low
K11
=
f low
K12

2
γz ,
m
1
γ2
=
ηz + z
m
σz

(3.28)
,

f low
K55
=
f low
K56

(3.29)
(3.30)

.

(3.31)

The K f orce terms (Eq. 3.26), which incorporate the internal Coulomb repulsion felt
by the electrons, are harder to solve due to the interaction potential Φ but, nevertheless,
substituting Eq. 3.15 and manipulating the resulting integrals allows them to be expressed in
terms of the pulse parameters. We will skip the details of this calculation as it is somewhat
cumbersome and tedious and limit ourselves to reporting the relevant results (for details
see [37]):

f orce
f orce
K11
= K55
= 0,

2γT f orce
K
,
σT 12
2γz f orce
=
K
,
σz 56
1
N e2
=
LT (ξ),
√
4π 0 6 σT π

f orce
K22
=

(3.33)

f orce
K66

(3.34)

f orce
K12

f orce
K56
=

where

0

(3.32)

1 N e2
Lz (ξ),
√
4π 0 6 σz π

(3.35)
(3.36)

is the vacuum permittivity, N is the number of electrons in the pulse, e and m are

53

the electron charge and rest mass, respectively, and

with ξ =

LT (ξ) =

ξ 2 L(ξ) − ξ
3
L(ξ) +
,
2
1 − ξ2

(3.37)

Lz (ξ) =

3ξ 2
[ξL(ξ) − 1] ,
ξ2 − 1

(3.38)

σz /σT > 0 and

L(ξ) =

1
2

π
0


√


arcsin 1−ξ 2

√

2

dθ
1−ξ
=
√
1 + ξ sin θ 
2


 ln(ξ+
√ ξ −1)
ξ 2 −1

for 0 ≤ ξ ≤ 1
(3.39)
for ξ > 1.

Note that for ξ = 1, L(1) = LT (1) = Lz (1) = 1. The functions L(ξ), LT (ξ) and Lz (ξ) (Fig.
3.2) describe the effect of the Coulomb interactions on the pulse’s time evolution and are
the only coupling in our model between the transverse and longitudinal degrees of freedom.

3
L(ξ)
LT(ξ)
Lz(ξ)

2.5
2
1.5
1
0.5
0
0

0.5

1

1.5

2

2.5

3

ξ = (σz/σT)1/2
Figure 3.2: Functions L(ξ) (Eq. 3.39), LT (ξ) (Eq. 3.37) and Lz (ξ) (Eq. 3.38) as a function
of ξ = σz /σT .

54

Once the Kijf low and Kijf orce terms are known, we can use Eqs. 3.20 and 3.21 to derive
differential equations for the parameters describing the pulse, σT,z , ηT,z and γT,z as

dσT (z)
∂
f low
= χ11(55) = K11(55)
,
(3.40)
dt
∂t
dγT (z)
∂
f low
f orce
= χ12(56) = K12(56)
+ K12(56)
,
(3.41)
dt
∂t
γT2 (z) ∂
γT2 (z) f low
dηT (z)
∂
2 ∂
2
f low
= χ22(66) −
χ12(56) + 2
χ11(55) = −
K12(56)
.
+ 2 K11(55)
dt
∂t
σT (z) ∂t
σT (z) ∂t
σT (z)
σT (z)
(3.42)
Substituting the previously calculated terms gives the final result:

dσi
2
= γi ,
dt
m
1
dγi
γ2
=
ηi + i
dt
m
σi

(3.43)
+

1 N e2
Li (ξ),
√
4π 0 6 σi π

dηi
2γi ηi
=−
,
dt
mσi

where i = T or z. Note that the quantity

(3.44)
(3.45)

√
σi ηi , corresponding to the pulse emittance, is

conserved in this formulation.
To check our calculation we consider the limit of zero Coulomb force, in which Eqs.
3.43-3.45 reduce to

dσi
2
= γi ,
dt
m
dγi
1
γ2
=
ηi + i
dt
m
σi
dηi
2γi ηi
=−
.
dt
mσi

(3.46)
,

(3.47)
(3.48)

This system of equations can be easily solved by noting that d2 γ/dt2 = 0 and by integrating,
55

with the result

ηi (0) 2
t + σi (0),
m
ηi (0)
γi (t) =
t,
m
σi (0)
ηi (t) = 2 2
.
t /m + σi (0)/ηi (0)

σi (t) =

(3.49)
(3.50)
(3.51)

These equations, together with Eq. 3.15, are an exact solution of Eq. 3.14, which in the
absence of the Coulomb force reads

∂f (r, p; t) p
∂f (r, p; t)
=−
· .
∂t
∂r
m

(3.52)

In the general case which includes the Coulomb force, no analytic solution of Eqs. 3.433.45 exists and one must numerically integrate the equations.

3.1.2

Lenses and RF cavities

So far we have described the time evolution of an electron pulse with no applied forces, but
for the Analytic Gaussian model to be useful in describing the UEM column, we need to
incorporate the effect of the optical elements (lenses and RF cavities) and we do so following
the work by Berger and Schroeder [39].
Rewriting the differential equation for the covariance matrix (Eq. 3.24) with an additional
external position dependent force term gives

∂
χi,j = Kijf low + Kijf orce + Kijext ,
∂t

56

(3.53)

where Kijf low and Kijf orce are defined in Eqs. 3.25 and 3.26 and

Kijext =

1
N

→
−
∂f (u; t)
du ui uj F ext (r) ·
,
∂p

(3.54)

→
−
where u = (x, px , y, py , z, pz ), r = (x, y, z), p = (px , py , pz ) and F ext (r) is the external force.
In this case Eqs. 3.43-3.45 become

dσT (z)
2
ext
= γT (z) + K11(55)
,
dt
m
γT2 (z)
dγT (z)
1
=
ηT (z) +
dt
m
σT (z)

(3.55)
+

1
N e2
ext
LT (z) (ξ) + K12(56)
,
√
4π 0 6 σT (z) π

2
2γT (z) ηT (z) γT (z) ext
dηT (z)
=−
+ 2 K11(55) ,
dt
mσT (z)
σT (z)

(3.56)
(3.57)

where we have made use of the fact that

ext
ext
K12(56)
= K21(56)
=

σT (z) ext
K
.
2γT (z) 22(66)

(3.58)

Using these equations, any optical element can be added into our system by specifying its
→
−
ext
ext
.
and K12(56)
force on the electron pulse, F ext (r), and by calculating K11(55)
In the transverse direction the electron pulse can be focused by using magnetic lenses.
Taking a perfect parabolic lens, the external force is linear with distance

→
−
F lens (r) = −MT xˆ
x − MT y yˆ,

(3.59)

where MT characterizes the strength of the lens and is positive for pulse focusing. Performing
ext
ext
ext
the integrals it is easy to see that K11(55)
= 0, K56
= 0 and K12
= MT σT , so that the

57

equations describing the system read

dσi
2
= γi ,
dt
m
dγi
1
γ2
=
ηi + i
dt
m
σi

(3.60)
+

1 N e2
Li (ξ) + Mi σi ,
√
4π 0 6 σi π

dηi
2γi ηi
=−
,
dt
mσi

(3.61)
(3.62)

with i = T or z and Mz = 0. To specify the physical region in which the lenses act, we take

Mi → Mi exp −

z − zlens
Llens /2

2n

,

(3.63)

where z is the position of the pulse center of mass in the lab frame, zlens is the position of the
center of the lens, Llens is the length of the lens and n is a super-Gaussian order parameter
that sets the sharpness of the lens edge, where n = 1 corresponds to a Gaussian profile and
n → ∞ is a top-hat.
To incorporate RF cavities into our treatment, we follow a similar procedure considering,
for simplicity, a cylindrical (pillbox) cavity operating in the TM010 mode. In this case the
electric and magnetic fields inside the cavity may be written in cylindrical coordinates and
to first order in r as [69]

Ez (r, z, z ) = E0 (z ) sin
Er (r, z, z ) = −

Ω(z − z)
+φ ,
v0

r ∂E0 (z )
sin
2 ∂z

r
Ω
Bφ (r, z, z ) = − E0 (z ) 2 cos
2
c

Ωz
+φ ,
v0
Ωz
+φ ,
v0

(3.64)
(3.65)
(3.66)

where Ω and φ are the RF field frequency and phase, so that the cavity will act as a pulse
58

compressor if φ

0, while for φ

±π/2 it will accelerate the pulse. z and v0 are the position

and velocity of the pulse center of mass in the cavity, which extends between z = −d/2 and
z = d/2, and r and z correspond to the position of an electron in the pulse with respect to
the center of mass. E0 (z ) is the field amplitude which we write in a form similar to that of
the lens strength as
E0 (z ) = E0 exp −

z − zrf
d/2

2n

,

(3.67)

where zrf is the position of the center of the cavity and n is the super-Gaussian order
parameter that sets the sharpness of the edge of the RF cavity field.
In the longitudinal direction, if the spatial dimension of the pulse is small compared to
the cavity size, Ωz << v0 , we can expand Eqs. 3.64 as

Ez (r, z, z ) = E0 (z ) sin

Ωz
Ωz
+φ −
cos
v0
v0

Ωz
+φ
v0

,

(3.68)

which shows that the RF cavity both accelerates the center of mass of the pulse (first term on
the RHS) and compresses (or expands) it (second term on the RHS). Since the accelerating
component is constant with respect to the pulse coordinates (r,z), it will not contribute
to the force but it must be included in the calculation of the center of mass velocity. The
resulting force on the pulse is therefore

Fz = −

eE0 (z )Ω
cos
v0

Ωz
+ φ z zˆ.
v0

(3.69)

In the transverse direction we substitute Eq. 3.67 into Eq. 3.65 and take the derivative,

59

with the result that

2n
Er (r, z, z ) = E0 (z )r
d

z − zrf
d/2

2n−1

sin

Ωz
+φ .
v0

(3.70)

The force in the radial direction can then be written as

Fr = eEr (r, z, z ) − ev0 Bφ (r, z, z )
= eE0 (z )

2n
d

z − zrf
d/2

(3.71)

2n−1

sin

Ωz
Ωv0
+ φ + 2 cos
v0
2c

Ωz
+φ
v0

rˆ
r.

(3.72)

ext
= 0,
From Eqs. 3.69 and 3.72, we can calculate the terms entering our equations as K11(55)
ext
= MTrf(z) σT (z) , with
and K12(56)

MTrf = eE0 (z )
Mzrf = −

3.1.3

2n
d

z − zrf
d/2

eE0 (z )Ω
cos
v0

2n−1

sin

Ωz
Ωv0
+ φ + 2 cos
v0
2c

Ωz
+φ .
v0

Ωz
+φ
v0

,

(3.73)
(3.74)

Photoemission gun

The last element we need to incorporate into our UEM column is the photoemission gun.
In the initial part of the microscope, inside the gun, the electrons are accelerated under a
constant electric field. Therefore particles in the tail of the pulse spend more time in the gun
region and are accelerated to a higher velocity. Overall this has the effect of decreasing the
chirp of the pulse as it exits the gun, and we need to incorporate this effect into our model.
We represent the field in the gun region with a step function, Egun = Θ(z − d) where d is
the spacing between cathode and anode and Θ is the Heaviside theta function. In the pulse
center of mass frame of reference, in which our Analytic Gaussian equations are formulated,
60

we have two cases: if the pulse center of mass is still inside the region with the electric
field, z < d, the portion of the pulse extending outside the anode (z > d − z ) will feel a
force F = −eEgun . On the other hand if z > d, the portion of the pulse before the anode,
z < d − z will be subject to a force F = eEgun . This is summarized in Fig. 3.3 and Eq.
3.75.
z’ < d

z’ > d

(a)

(b)
-eEgun

z’

d

eEgun

d

z

z’

z

Figure 3.3: Effect of stepwise electric field of photoemission gun: (a) z < d, (b) z > d





Fz = −eEgun zˆ for z < d,

z >d−z




Fz = eEgun zˆ

z < z − d.

(3.75)

for z > d,

Following the procedure used to incorporate the lenses and RF cavity, we write Eq. 3.54 as

Kijgun =

1
N

du ui uj Fz

∂f (u; t)
.
∂pz

(3.76)

Inserting the expression for f (u; t) (Eq. 3.15) and performing the derivative, gives

Kijgun =

1
N

du ui uj Fz

1
ηz

pz −

γz
z f (u; t),
σz

(3.77)

with u = (x, px , y, py , z, pz ). In the previous section we showed how any external force could
ext
ext
ext
ext
be added by simply calculating the elements K11
, K12
, K55
and K56
. In this case it is

61

gun
gun
gun
vanish, so the only term of interest is
and K12
, K55
clear that the terms K11

gun
K56
=

=

1
N

du z pz Fz

1
2π σz ηz3

=√

1
2πσz

1
ηz

pz −

γz
z f (u; t)
σz

dzdpz z pz Fz pz −
dz Fz z exp −

z2
2σz

(3.78)

z2
[pz − (γz /σz )z]2
γz
z exp −
−
σz
2σz
2ηz

,

(3.79)
(3.80)

where the second line was obtained after integration over x, y, px , py and in the third line we
integrated over pz as well. Now we have two possibilities, if z < d, then the force is nonzero
only for z > d − z and

gun
K56
= −eEgun √

= −eEgun

1
2πσz

∞

dz z exp −
d−z

σz
(d − z )2
exp −
2π
2σz

z2
2σz

.

(3.81)
(3.82)

A calculation for z > d and z < z − d gives the same final result. We conclude that exiting
the gun region has the effect of a reduction in the chirp γz of the pulse proportional to

√

σz

multiplied by a spatial profile that decays as exp[−(d − z )2 /(2σz )].

3.1.4

Relativistic transformations

The equations we have written so far describe the electron pulse in its center of mass frame
of reference, while our observations are performed in the laboratory frame of reference. If
the center of mass of the pulse is moving at non-relativistic speeds (v0 /c

1, where v0 is the

magnitude of the COM velocity and c is the speed of light), a simple classical transformation

62

can be used:
flab (r, p; t) = f (r − r0 − v0 t, p − mv0 ; t),

(3.83)

where r0 and v0 are, respectively, the initial center of mass position and velocity.
The addition of an external constant electric field E can also be handled in a simple way,
since the transformation to the lab frame of the distribution f (r, p; t) simply becomes

flab (r, p; t) = f (r − r0 − v0 t, p − mv0 − eEt; t).

(3.84)

In the context of typical UEM systems, where the pulse is initially accelerated with
extraction fields ranging from 0.1 MV/m to 100 MV/m, resulting in velocities of the center
of mass in the z-direction (taken parallel to the column axis) that are some fraction of the
speed of light, these transformations are not applicable as it is necessary to take relativistic
effects into account. In the center of mass frame the electrons are not moving relativistically
with respect to each other, so the only modification to the Analytic Gaussian model comes
from the transformations that connect the COM frame to the lab frame of reference.
In general, the transformations from a coordinate frame at rest (F) to a frame moving at
a velocity v (F’) can be easily written if one first defines the rapidity β [70, 71] as a function
of the velocity v:
v
= tanh β,
c

(3.85)

where we assume v to be parallel to the z-axis. In this case the Lorentz transformations are

63

simply written as

z = z cosh β + ct sinh β,

(3.86)

ct = ct cosh β + z sinh β,

(3.87)

where z, t are coordinates in the rest frame (F), while the primed quantities denote variables
in the moving frame (F’). Note that in this notation the Lorentz factor γL = cosh β. For a
generic particle moving at a speed u = c tanh α with respect to the rest frame F, both energy
E and momentum p can be expressed as

E = mc2 cosh α,

(3.88)

p = mc sinh α,

(3.89)

and the transformation of these quantities to the frame F’ moving at a speed v is conveniently
handled due to the additivity of the rapidity, resulting in

v = c tanh(α + β),

(3.90)

E = mc2 cosh(α + β),

(3.91)

p = mc sinh(α + β).

(3.92)

Now, assuming the electron pulse in our model has been accelerated to a velocity v0 =
(0, 0, v0 ), where v0 /c > 0, and that no external fields are acting on it, the transformation
from the distribution in the COM frame to the lab frame of reference is conveniently given

64

by [70]
flab (x, y, z, px , py , pz ; t) = f (x , y , z , px , py , pz ; t ),

(3.93)

with

v0
= tanh β,
c
x = x,
y = y,
z = z cosh β − ct sinh β,
px = p x ,
py = p y ,
pz = mc sinh(α − β),
ct = ct cosh β − z sinh β,

where α is defined by pz = mc sinh α.
On the other hand for a constant electric field, E0 = (0, 0, Ez ), the COM frame accelerates
with respect to the lab frame with an acceleration in the z-direction equal to a = eEz /m
and the transformation between the two coordinate systems now reads

z + z0 = (z + z0 ) cosh [β(t )] ,

(3.94)

ct = (z + z0 ) sinh [β(t )] ,

(3.95)

where, again, z, t are coordinates in the lab frame while the primed quantities denote vari-

65

ables in the pulse COM frame and we have defined
c2
,
a
at
β(t ) =
.
c
z0 =

(3.96)
(3.97)

In the limit of zero acceleration, these equations reduce to the transformation previously
written.
Using these equations, it is possible to compute the time evolution of the pulse in its COM
frame of reference using the Analytic Gaussian equations and then transform the resulting
distribution into the lab frame of reference. More importantly, as we are interested in the
pulse parameters σi , γi , and ηi , with i = T or z, we know that they are simply transformed
as

σTlab = σT ,

(3.98)

σzlab = σz / cosh β = σz /γL2 ,

(3.99)

γTlab = γT ,

(3.100)

γzlab = γz ,

(3.101)

ηTlab = ηT ,

(3.102)

ηzlab = ηz cosh β = ηz γL2 ,

(3.103)

where β is the rapidity of the COM with respect to the lab frame, γL is the Lorentz factor
and the equations we wrote simply represent the well known length contraction and time
dilation phenomena in relativistic mechanics (note that in our notation σi has units of length
squared and, similarly, ηi is in units of momentum squared).
66

In addition to this, it is also necessary to consider the relativistic transformations of the
applied electric and magnetic fields. Taking the fields in the lab frame to be known, E and
B, their transformation to the pulse frame is given by

E⊥ = γL (E⊥ + v0 × B),

(3.104)

Ez = Ez ,

(3.105)

B⊥ = γL (B⊥ −
Bz = Bz ,

1
v0 × E),
c2

(3.106)
(3.107)

where γL is the Lorentz factor, v0 = (0, 0, v0 ) is the velocity of the pulse center of mass, E
and B are the electric and magnetic fields in the lab frame, and the primed quantities denote
variables in the pulse frame.
In the COM frame the average velocity is equal to zero, so the force acting on the pulse
center of mass is given by F = eE . The longitudinal electric field is not modified so we
have that Fz = eEz = eEz = Fz . To consider the effect in the transverse direction, it is
convenient to think in terms of cylindrical coordinates. For a magnetic lens, Bφ = 0, E = 0,
so that Er = −γL v0 Bφ and Fr = −qγL v0 Bφ = γL Fr , where Fr = −qv0 Bφ is the force in the
lab frame. In our formalism this corresponds to an increase of the lens strength equal to γL .
For an RF cavity, Er = 0, Bφ = 0 and Er = γL (Er − v0 Bφ ), resulting in a similar effect.
These contributions can be incorporated in our formulation by renormalizing the transverse

67

lens and RF cavity strengths in the following way:

MTlens → γL MTlens ,

(3.108)

Mzlens = 0,

(3.109)

MTrf → γL MTrf ,

(3.110)

Mzrf = Mzrf .

(3.111)

The field in the photoemission gun, which is entirely in the z direction, is not affected as
fields parallel to the relativistic motion are not transformed. However, our derivation of the
effect of the gun region resulted in a term ∝ exp [−(dgun − z )2 /(2σi )], where dgun and z are
the length of the photoemission gun and the position of the pulse center of mass in the lab
frame. It becomes clear that the width of the pulse appearing in this expression should also
be taken in the lab frame, meaning that we perform the substitution σi → σi /γL2 .
For reference we now report the equations governing the Analytic Gaussian model with
the inclusion of the optical elements discussed in the text and the effects of relativistic
motion. Note that the units for the parameters are: σi units of length squared, ηi units of
momentum squared and γi units of action (length x momentum).

dσi
2
= γi ,
dt
m
dγi
1
γ2
=
ηi + i
dt
m
σi

(3.112)
+

√
1 N e2
(dgun − z )2
Li (ξ) + (Milens + Mirf )σi − Migun σi exp −
,
√
4π 0 6 σi π
2σi /γL
(3.113)

dηi
2γi ηi
=−
,
dt
mσi

(3.114)

68

with i = T or z and

ξ=

σz /σT ,

LT (ξ) =

(3.115)

3
ξ 2 L(ξ) − ξ
L(ξ) +
,
2
1 − ξ2

(3.116)

3ξ 2
Lz (ξ) = 2
[ξL(ξ) − 1] ,
ξ −1

√


arcsin 1−ξ 2

√

for 0 ≤ ξ ≤ 1
1−ξ 2
L(ξ) =
√


ln(ξ+ ξ 2 −1)

 √
for ξ > 1,

(3.117)

(3.118)

ξ 2 −1

MTlens = γL MT exp −

2n

z − zlens
Llens /2

,

(3.119)

Mzlens = 0,
MTrf

(3.120)
z − zrf
d/2

2n
= γL eE0 (z )
d

Mzrf = −

eE0 (z )Ω
cos
v0

E0 (z ) = E0 exp −

2n−1

sin

Ωz
Ωv0
+ φ + 2 cos
v0
2c

Ωz
+φ ,
v0

z − zrf
d/2

Ωz
+φ
v0

,

(3.121)
(3.122)

2n

,

(3.123)

MTgun = 0,

(3.124)

1
Mzgun = eEgun √ .
2π

(3.125)

See previous text for further details of the notation used.

69

3.2

Methods

The set of equations we have derived (Eqs. 3.112-3.125) describing all the optical elements
in our column, have to be numerically integrated to give the time evolution of the electron
pulse inside the microscope. To do so we developed a program written in Fortran 90 that
performs this numerical integration using a fourth order Runge Kutta method with adaptive
step size control, implemented following the procedure outlined in [72]. Given the positions
and properties of the elements one wishes to use, it is then possible to calculate the beam
size, chirp and momentum spread at the target position.
In order to provide information to guide the development of the electron microscope, an
optimization routine was also implemented to adjust the strengths of the magnetic lenses
and RF cavity in order to focus both the longitudinal and transverse pulse components at the
target and minimize the pulse size at the focal point. This was done using the Levenberg–
Marquardt algorithm [72], a quasi-Newton optimization method which allows to minimize a
given objective function F that depends on a set of parameters p. In our case the components
of p correspond to the RF cavity and magnetic lens strengths and what we want to achieve
is a focal spot in which both the transverse and longitudinal widths have a minimum value.
To do this we choose to use the following objective function

F = c1

zmin,z − zmin,T
zmin,z

2

+ c2

σmin,z
σ0,z

+ c3

σmin,T
σ0,T

,

(3.126)

where c1 , c2 , c3 are the relative weights of each term, σmin,z(T ) is the minimum pulse width
squared (following the notation used in our model) and zmin,z(T ) is the corresponding position
in the longitudinal (transverse) direction. σ0,z , and σ0,T are the initial pulse widths.
Given our function F (p) and a set of parameters pi , the Levenberg–Marquardt algorithm
70

updates these parameters to minimize F (p) in the following way

pi+1 = pi − (H + λHdiag )−1 d,

(3.127)

where

Hi,j =

∂F (p) ∂F (p)
,
∂pi
∂pj

d = F (p)

∂F (p)
.
∂p

(3.128)
(3.129)

H is an approximation of the true Hessian matrix obtained by multiplying the first order
derivatives and Hdiag is a matrix formed with the diagonal elements of H. The parameter
λ, which is adjusted throughout the optimization procedure, determines the behavior of the
algorithm: for large λ the above equation approaches the steepest descent method, while
for small λ it approaches the Gauss-Newton algorithm. In practice we adjust λ so that far
from the minimum we use the former, while as we get closer to a minimum of F (p), we
incorporate more of the latter. We stop the minimization routine after performing a certain
number of steps (typically we choose N=3) where the fractional change in F (p) is less than
10−3 . Once convergence has been reached, an evaluation of the matrix H gives information
on the covariance of the optimization parameters.

3.3

Results

In this section we present the results obtained using the Analytic Gaussian model and the
code we developed. We start by looking at the time evolution of an electron pulse in a
drift region and compare the Analytic Gaussian results to N particle simulations done using
71

the COSY-MLFMM code described in chapter 2. We then proceed to analyze the effect of
each optical element on the electron pulse and then combine them to simulate the whole
microscope column.

3.3.1

Drift

It is clear that our simple mean field Analytic Gaussian model cannot capture the details
of the dynamics of an electron pulse, since it assumes a Gaussian pulse shape at all times
and neglects any nonlinear contributions of the optical elements. Nevertheless, we expect it
to be reasonably accurate in those sections of the column where the higher order electron–
electron correlations and aberrations can be ignored. In particular, since there are long drift
segments between the optical components of the microscope, we want to understand how well
our model performs in describing them. Having observed in our photoemission simulations
that there are two regimes of interest, below and above the virtual cathode (VC) limit, we
compare the predictions of the Analytic Gaussian model in both of these cases by fitting the
pulse profile at the end of the photoemission process (t =120 ps) with the Gaussian profile
in Eq. 3.15 and then calculating the time evolution over a period of 50 ps using both codes.
Below the VC and for a Gaussian laser profile, even though the photoemitted pulse
displays some nonlinear distortions both in the transverse and longitudinal directions, the
Analytic Gaussian model is able to describe the evolution of the pulse phase space (Fig.
3.4) and capture the time evolution of the pulse parameters remarkably well (Fig. 3.5),
with fractional errors less than 1.5 %. The pulse momentum spread is calculated from the
Analytic Gaussian parameters as
σpi = ηi +

72

γi2
,
σi

(3.130)

where i = T or z and ηi .
(a)

2

(b)

2
t=50 ps

1

px (10-3 MeV/c)

px (10-3 MeV/c)

t=0.0 ps

0

-1

-2

1

0

-1

-2
-300 -200 -100 0 100 200 300
x (µm)

3

pz (10-3 MeV/c)

2

(d)
t=0.0 ps

3
2

1

pz (10-3 MeV/c)

(c)

-300 -200 -100 0 100 200 300
x (µm)

0
-1
-2
-3

t=50 ps

1
0
-1
-2
-3

-4

-4
-100

0
z (µm)

100

-200

-100

0
z (µm)

100

200

Figure 3.4: Color map of the charge distribution of an electron bunch below the virtual
cathode limit (Fa =1MV/m, Ne =106 , Gaussian laser profile) in the rest frame of the bunch
projected onto the px -x plane [(a), (b)] and pz -z plane [(c),(d)]. Also shown is the time
evolution of the Analytic Gaussian model for the pulse, with the black line indicating the
contour where the intensity drops to 0.5 times its peak value.

Since our photoemission simulations showed that it is advantageous to use an Elliptical transverse laser profile, in Figs. 3.6 and 3.7 we compare the N particle and Analytic Gaussian simulations with the electron pulse generated using an elliptical laser profile
[I(r) ∝

1 − (r/R)2 ]. Our Analytic Gaussian model is able to accurately simulate the

time evolution in this case as well, and in fact, due to the more linear nature of the phase

73

(b) 140
√σT (µm)


1
0.75
0.5

130

70
Longit. (z)
AG
N-part.

120

0.25

50

110
20 30
t (ps)

40

50

(c) 12

0
(d)

11
10

11
10

9

9
Longit. (z)
AG
N-part.

8
7

γz (10-2 µm MeV/c)

12
Transv. (T)
AG
N-part.

8
7

6

6
0

10

20 30
t (ps)

40

10

20 30
t (ps)

40

50

9

15
Transv. (T)
AG
N-part.

8

14

Longit. (z)
AG
N-part.

7

T

10

√σp (10-4 MeV/c)


0
γT (10-2 µm MeV/c)

60

13

6

50

z

0

80

√σp (10-4 MeV/c)


z (mm)

90
Transv. (T)
AG
N-part.

AG
N-part.

√σz (µm)


(a)

12
0

10

20 30
t (ps)

40

50

Figure 3.5: Time evolution of pulse parameters below the virtual cathode limit (Fa =1MV/m,
Ne =106 , Gaussian laser profile) comparing N-particle simulations and the Analytic Gaussian
model. Panel (a) shows the pulse center of mass position, panels (b)-(d) show the transverse
√
(T) and longitudinal (z) pulse width σi [panel (b)], chirp γi [panel (c)] and momentum
√
spread σpi [panel (d)].
space, the fractional errors in the parameters are smaller than for the previous case, with all
parameters presenting a difference of less than 1.0 %.
On the other hand, above the virtual cathode limit, the pulse is very strongly nonlinear
(Fig. 3.8) and it is clear that a Gaussian fit of this pulse shape does not reflect its complex
dynamics. The time evolution of the pulse parameters, using a Gaussian laser profile, is
shown in Fig. 3.9. While the spatial width of the pulse is still reasonably captured, with a
fractional error that after 50ps is less than 1%, the momentum distribution is not accurately
described. In the transverse direction the fractional errors for γT and

√

σpT are between

3-4%, while in the longitudinal direction the error is 5-6% at t=50ps. As a result of these,
for longer time simulations, we expect that also the spatial width σi would deviate from the
N-particle prediction. Similar results hold for the Elliptical laser profile.

74

(a)

2

(b)

2
t=50 ps

1

px (10-3 MeV/c)

px (10-3 MeV/c)

t=0.0 ps

0

-1

-2

1

0

-1

-2
-300 -200 -100 0 100 200 300
x (µm)

3

pz (10-3 MeV/c)

2

(d)
t=0.0 ps

3
2

1

pz (10-3 MeV/c)

(c)

-300 -200 -100 0 100 200 300
x (µm)

0
-1
-2
-3

t=50 ps

1
0
-1
-2
-3

-4

-4
-100

0
z (µm)

100

-200

-100

0
z (µm)

100

200

Figure 3.6: Color map of the charge distribution of an electron bunch below the virtual
cathode limit (Fa =1MV/m, Ne =106 , Elliptical laser profile) in the rest frame of the bunch
projected onto the px -x plane [(a), (b)] and pz -z plane [(c),(d)]. Also shown is the time
evolution of the Analytic Gaussian model for the pulse, with the black line indicating the
contour where the intensity drops to 0.5 times its peak value.
From these results we conclude that below the virtual cathode limit, the Analytic Gaussian is indeed able to describe the evolution of the pulse with a high degree of fidelity even
though the N-particle distribution is not truly Gaussian. Above the VC limit, on the other
hand, this approximation breaks down and major differences start showing already after
50ps.

75

(b) 140
√σT (µm)


1
0.75
0.5

130

70
Longit. (z)
AG
N-part.

120

0.25
0

80

60
50

110
50
12
11
10
9
8
7
6
5

Transv. (T)
AG
N-part.
Longit. (z)
AG
N-part.
0

10

20 30
t (ps)

40

0

10

20 30
t (ps)

40

50

(d)

14
Transv. (T)
AG
N-part.

8

7

13

Longit. (z)
AG
N-part.

z

40

T

20 30
t (ps)

√σp (10-4 MeV/c)


γT (10-2 µm MeV/c)

(c) 12
11
10
9
8
7
6
5

10

γz (10-2 µm MeV/c)

0

√σp (10-4 MeV/c)


z (mm)

90
Transv. (T)
AG
N-part.

AG
N-part.

√σz (µm)


(a)1.25

12

6

50

0

10

20 30
t (ps)

40

50

Figure 3.7: Time evolution of pulse parameters below the virtual cathode limit (Fa =1MV/m,
Ne =106 , Elliptical laser profile) comparing N-particle simulations and the Analytic Gaussian
model. Panel (a) shows the pulse center of mass position, panels (b)-(d) show the transverse
√
(T) and longitudinal (z) pulse width σi [panel (b)], chirp γi [panel (c)] and momentum
√
spread σpi [panel (d)].

3.3.2

Photoemission gun

Next we turn our attention to the various elements that compose the microscope, starting
from the effect that crossing the anode plane of the photoemission gun, corresponding to
the end of the extraction field, has on the pulse properties. We consider varying numbers of
electrons below the virtual cathode limit, which from our previous photoemission simulations
we estimate to be about Ne = 1.5 · 107 for Fa =1MV/m. As the transverse components are
not directly affected by this so-called acceleration gap, we focus on the longitudinal pulse
parameters in Fig. 3.10. The slower electrons in the tail of the pulse spend more time in
the acceleration region and are accelerated more, resulting in a reduction of the pulse chirp,
as seen in the Figure. This effect is seen to be more step-like for lower number of electrons
(Ne = 105 ) due to the smaller pulse size σz compared to the Ne = 107 case.
76

(a)

5

(b)

5
t=50 ps

2.5

px (10-3 MeV/c)

px (10-3 MeV/c)

t=0.0 ps

0

-2.5

-5
-600

-300

0
x (µm)

300

2.5

0

-2.5

-5
-600

600

(c)

0
x (µm)

300

600

(d)
t=0.0 ps

t=50 ps
10
pz (10-3 MeV/c)

10
pz (10-3 MeV/c)

-300

0

-10

-20

0

-10

-20
-600 -300

0
300
z (µm)

600

-1000 -500

0
500
z (µm)

1000

Figure 3.8: Color map of the charge distribution of an electron bunch above the virtual
cathode limit (Fa =1MV/m, Ne =108 , Gaussian laser profile) in the rest frame of the bunch
projected onto the px -x plane [(a), (b)] and pz -z plane [(c),(d)]. Also shown is the time
evolution of the Analytic Gaussian model for the pulse, with the black line indicating the
contour where the intensity drops to 0.5 times its peak value.

3.3.3

Lens

We now analyze the treatment of magnetic lenses with the Analytic Gaussian model. Each
lens is described by three parameters: the strength M0 , the spatial width Llens , and an order
parameter n (n = 1 corresponds to a Gaussian profile). Given the magnetic field data for
the lenses used in the microscope, we use Eq. 3.119 with γL = 1 to perform a fit and obtain
the following values: n = 1 and Llens = 2.157mm. M0 , proportional to the current in the

77

0.75
0.5
0.25
0

250

500
Longit. (z)
AG
N-part.

225

40

50

400

1

8
Transv. (T)
AG
N-part.

0.75

6

0.5
Longit. (z)
AG
N-part.

0.25

300
0

4

10

20 30
t (ps)

40

50

(d)

105
Transv. (T)
AG
N-part.

26

24

100

Longit. (z)
AG
N-part.

T

20 30
t (ps)

√σp (10-4 MeV/c)


10

γz (µm MeV/c)

γT (µm MeV/c)

600

200
0

(c)

275

95

22
0

10

20 30
t (ps)

40

50

z

√σT (µm)


z (mm)

1

700
Transv. (T)
AG
N-part.

√σz (µm)


(b) 300
AG
N-part.

√σp (10-4 MeV/c)


(a)

90
0

10

20 30
t (ps)

40

50

Figure 3.9: Time evolution of pulse parameters above the virtual cathode limit (Fa =1MV/m,
Ne =108 , Gaussian laser profile) comparing N-particle simulations and the Analytic Gaussian
model. Panel (a) shows the pulse center of mass position, panels (b)-(d) show the transverse
√
(T) and longitudinal (z) pulse width σi [panel (b)], chirp γi [panel (c)] and momentum
√
spread σpi [panel (d)].
magnetic lens, is left as a variable to adjust the focal length, as shown in Fig. 3.11 where
we plot the variation of the transverse pulse parameters due to a lens at z=45mm. In our
derivation, the lens affects dγT /dt directly with a term proportional to M0 σT , which is seen
in the change of γT in Fig. 3.11(c). The reduction in γT results in a focusing of σT due to
the coupling between the parameters and, since the emittance =

√

σi ηi , is conserved in our

model, this results in a defocusing of ηT . As the strength of the lens is decreased, the focal
distance increases and so does the pulse size at the focal point.
For a given lens, for which we take M0 = 5, the focusing power depends on the number of
electrons and the emittance of the pulse. Given the same initial conditions [Fig. 3.12(a)], the
higher the number of electrons, the harder it is to focus the pulse. In our data this is seen in
the increased focal length and spot size at the focal point in the case of the Ne = 107 pulse.

78

(b) 104

(a) 0.5

√σz (µm)


v0/c

0.4
0.3
0.2

Ne=105
Ne=1067
Ne=10

0.1

103
102

Ne=105
Ne=1067
Ne=10

0
0

30

60

90

0

30

z (mm)

90

(d) 100

(c)
5

5

√ηz (me µm/ps)


Ne=106
Ne=107
Ne=10

20
γz/γ0

60
z (mm)

10

0
0

30

60

90

z (mm)

Ne=106
Ne=107
Ne=10

10-1
10-2
10-3

0

30

60

90

z (mm)

Figure 3.10: Longitudinal pulse parameters as a function of distance along the microscope
column for Fa =1MV/m, as the pulse crosses the anode plane (at z=45mm, indicated with
dashed line). Panel (a) shows the center of mass velocity as a fraction of the speed of light c,
√
while panels (b)-(d) show the variation of the longitudinal (z) pulse width σz [panel (b)],
√
chirp γz /γ0 , where γ0 is the initial chirp [panel (c)] and local momentum variance ηz [panel
(d)].
On the other hand, from our study of the photoemission process, we know that increasing
the number of electrons emitted results also in an increase in the pulse emittance. Therefore,
using the data from our photoemission simulations as the initial condition corresponds to
different initial values of the pulse emittance as Ne is varied. In our model pulses with a
higher emittance are focused more easily and to a smaller spot size than those with a lower
emittance, since focusing can be achieved at the expense of increasing the momentum spread
ηi . This results in pulses with a higher Ne , which have a higher emittance, being focused
better that those with lower Ne , as seen in Fig. 3.12(b), which shows that the focal length
of the lens is approximately the same for all the pulses considered with the Ne = 107 pulse

79

(a) 1

(b)

M0= 5
M0=10
M0=20

4

10
√σT (µm)


0.75
M(z)/M0

5

10

0.5
0.25

103
102
101

0
40 42.5 45 47.5 50
z (mm)
(c)

0

50

20

(d)

100
150
z (mm)
1

√ηT (me µm/ps)


γT/γ0

0
M0=5
M0=10
M0=20

-10
0

50

100 150
z (mm)

200

250

250

M0=5
M0=10
M0=20

10
10

200

100
10-1
10-2
10-3

0

50

100 150
z (mm)

200

250

Figure 3.11: (a) Fractional lens strength M (z)/M0 as a function of position. (b)-(d) Transverse pulse parameters as a function of distance along the microscope column for Ne =106 ,
with a magnetic lens of varying strength M0 at z=45mm (indicated with a dashed line).
√
Shown are the transverse (T) pulse width σT [panel (b)], chirp γz /γ0 [panel (c)] and local
√
momentum variance ηz [panel (d)].
displaying a smaller size at the focal spot than the other two.
These results lead to two observations: first they highlight the importance of correctly
simulating the photoemission process. The previous works using the Analytic Gaussian
model [37–39] did not consider the details of how the electrons are emitted and assumed
that the initial conditions would be the same regardless of the number of electrons in the
pulse. We have shown here how this is not the case and that varying the initial emittance
changes the focusing obtained in this model.
Our second consideration is about the limits of the Analytic Gaussian model, since so
far we have assumed that the lenses used are linear regardless of focusing strength and pulse

80

Fixed IC

N-part. IC

(a)

(b)
1000
√σT (µm)


√σT (µm)


1000

100
Ne=105
Ne=1067
Ne=10
10

100
Ne=105
Ne=1067
Ne=10
10

0

50

100 150 200 250 300
z (mm)

0

50

100 150 200 250 300
z (mm)

√
Figure 3.12: Transverse (T) pulse width σT as a function of distance along the microscope
column with a magnetic lens (M0 = 5) at z=45mm (indicated with a dashed line) for varying
numbers of electrons in the pulse. (a) Effect of lens given the same initial conditions. (b)
Effect of lens using the initial conditions from the photoemission simulations.
size. This approximation is only valid in the limit where the pulse size is much smaller
than the lens diameter. Checking the magnetic field calculated for the actual lens in the
microscope, we see that it retains a mostly Gaussian shape for r
for values of

√

2σT greater than this, or

√

0.8mm, meaning that

σT > 0.57mm, this approximation is not correct.

Comparing this value to the data shown in Fig. 3.12 shows that the width of the Ne = 107
pulse is well above this value before the lens, so that in a real system it would not be
focused to such a tight spot as our model predicts. To improve this, in principle, one could
include spherical aberrations in the Analytic Gaussian model by adding nonlinear external
forces, but it would disregard the mean field approach inherent in its derivation and cause
pulse breakup. In practice this means that while the AG model can give an estimate of the
properties of the column, these should be seen as best-case estimates and one should use
other simulation tools to capture the details of the dynamics, specially as the number of
electrons in each pulse is increased.

81

3.3.4

RF cavity

The last element to analyze in our column is the RF cavity, which is used to focus the pulse
in the longitudinal direction at the expense of a defocusing of the transverse components,
as seen in Figure 3.13. The parameters used for the RF cavity in our simulations are n =2,
Lrf =22.5mm, Ω = 2π· 1.013GHz and Er f =1.0-2.4MV/m, matching the parameters of the
cavity used in the actual microscope. In order to focus the pulse, the time dependent field
in the RF cavity needs to be zero when the center of mass of the pulse is at the center of the
cavity, so that the tail of the pulse is accelerated and the front is decelerated, with no net
change of the COM velocity [see Fig. 3.13(a)]. To do this it is necessary to match the phase
of the RF cavity with the time it takes the electron pulse to reach the center of the cavity,
which is not trivial to do in practice in a real microscope. Also, in addition to this, the time
of flight of the pulse through the cavity should be at most half of the oscillation period, to
prevent unwanted defocusing due to the oscillating field in the cavity. This sets a limit on
the maximum longitudinal pulse size that a cavity at a given frequency can compress, which
for our cavity and 100keV electrons can be estimated to be

√
σz

20mm.

As seen in the discussion on the magnetic lenses, the increased Coulomb repulsion present
in pulses with higher numbers of electrons, makes them harder to focus and this is true also
in the longitudinal direction (Fig. 3.14), resulting in a larger spot size at the focal point for
increasing Ne .

3.3.5

Column

So far we have discussed the optical elements individually, now we combine them to model
the proposed design for the Ultrafast Electron Microscope under development in our group.

82

0.5
v0/c
Erf/E0
Erf (MV/m)
2.4
1.0

103

0.25
0

0.38

√σi (µm)


v0/c

0.4

(b)

Erf/E0

(a)
0.42

-0.25

0.36

Erf (MV/m)
2.4
1.0
Long.
Transv.

101
100

-0.5
30

102

60

0

30

z (mm)
(c)

(d)

γi/γ0

20

√ηi (me µm/ps)


Erf (MV/m) Long.
Transv.
2.4
1.0

30

10
0
-10
0

30

60
90
z (mm)

120

1

10

Erf (MV/m)
2.4
1.0

100

120

150

Long.
Transv.

10-1
10-2
10-3

150

60
90
z (mm)

0

30

60
90
z (mm)

120

150

Figure 3.13: Pulse parameters as a function of distance along the microscope column for
Ne =106 with an RF cavity centered at z=45mm (indicated by the vertical dashed line) for
different RF field strengths. (a) Center of mass velocity as a fraction of speed of light c and
electric field as a fraction of the peak value E0 . (b)-(d) Transverse (T) and longitudinal (z)
√
√
pulse width σi [panel (b)], chirp γi /γ0 [panel (c)] and local momentum variance ηi [panel
(d)]. The continuous (dashed) lines indicate the longitudinal (transverse) components and
the colors correspond, respectively, to Erf = 2.4M V /m (black) and Erf = 1.0M V /m (red).
The positions and properties of the photoemission gun, magnetic lenses and RF cavity used
in the microscope are summarized in Table 3.1. The strengths of the lenses and of the electric
field in the RF cavity were used as variables in the optimization routine to control the focal
point position and spot size. The initial pulse parameters (number of electrons, center of
mass velocity and Analytic Gaussian parameters) were extracted from our simulations of the
photoemission process discussed in the previous chapter.
Results for the column optimized for Ne = 106 electrons generated with a Gaussian laser
pulse are shown in Fig. 3.15, where we plot the longitudinal and transverse pulse width

83

√

σi ,

(b) 104

(a)

103

103

Ne=105
Ne=1067
Ne=10
0

50

100
z (mm)

102

√σz (µm)


√σT (µm)


104

101
Ne=105
Ne=1067
Ne=10

100
10-1

150

0

50

100
z (mm)

150

√
√
Figure 3.14: (a) Transverse pulse width σT and (b) longitudinal pulse width σz as a
function of distance along the microscope column for varying Ne , with an RF cavity centered
at z=45mm (indicated by the vertical dashed line) with Erf =2.4MV/m.
Photoemission gun
Extraction field
Cathode position
Anode position
Magnetic lenses
Lens 1
Lens 2
Lens 3
RF cavity
Ω/ 2π = 1.013GHz
Sample

2.2MV/m
z=0mm
z=45mm
n=1, Llens =2.157mm
z=70.2mm
z=360.0mm
z=1066.0mm
n=2, Lrf =22.5mm
z=640.8mm
z=1080.5 ± 25.4 mm

Table 3.1: Parameters used in the simulation of the UEM system using the Analytic Gaussian
model.
chirp γi and local momentum variance

√

ηi as a function of position along the microscope.

The focal point is at 1088.62mm, with pulse size

√
√
σT = 0.92µm and σz = 2.27µm. Due

to the strength of the lens used to focus the pulse, for positions close to the focal point, the
transverse width diverges rapidly so accurate positioning of the sample at the focal spot is
crucial for good spatial resolution. The divergence of the longitudinal component is much
less severe.
We perform a similar optimization varying the number of electrons in the pulse and the
84

4

(a) 10

Long.
Transv.

3

√σi (µm)


10

2

10

101
100

γi (103 me µm2/ps)

(b)

√ηi (me µm/ps)


(c)

Long.
Transv.

1

0

-1
101

Long.
Transv.

100
10-1
10-2
10-3

0

200

400

600
z (mm)

800

1000

1200

Figure 3.15: Longitudinal (black) and transverse (dashed red) pulse parameters as a function of distance along the microscope column for Ne =106 . Shown are the (a) pulse width
√
√
σi , (b) chirp γi and (c) local momentum variance ηi . The region corresponding to the
photoemission gun is shown in blue, while that corresponding to the RF cavity is in green.
The dashed lines correspond to the magnetic lenses.
transverse profile of the photoemitting laser and summarize our results in Fig. 3.16 and
Table 3.2. We see that in our model the position of the focal spot is relatively insensitive to
the charge in the pulse. The longitudinal width

√
σz and the corresponding time resolution,

√
defined as δt = 2 2σz /v0 , are increased from 0.76µm (δt=13fs) for the 105 Gaussian electron
pulse to 17µm (δt=288fs) for the Gaussian pulse with 107 electrons. The transverse pulse
sizes are comparable for all the simulations reported using a Gaussian laser profile but, as

85

discussed in section 3.3.3, this is a result of neglecting aberrations which would certainly
play a role as the width of the pulse is increased. Using an elliptical laser profile to generate
the electron pulse offers some advantages in terms of a decreased transverse emittance,
which leads to a focal spot size that is half of what is predicted for the Gaussian case. In
the longitudinal direction the laser profile plays a less significant role, and the temporal
resolution is comparable to that obtained with a Gaussian laser.
5

6

(b) Gaussian - Ne=10

104

104

103

103
√σi (µm)


√σi (µm)


(a) Gaussian - Ne=10

102
101
Long.
Transv.

100
10-1

1050

1100
z (mm)

102
101

10-1

1150

104

104

103

103

101
Long.
Transv.

100
10-1

1050

1100
z (mm)

1050

1100
z (mm)

1150

(d) Elliptical - Ne=106

√σi (µm)


√σi (µm)


(c) Gaussian - Ne=107

102

Long.
Transv.

100

102
101
Long.
Transv.

100
10-1

1150

1050

1100
z (mm)

1150

Figure 3.16: Detail of the transverse and longitudinal pulse widths in the optimized column
around the focal spot for a Gaussian laser pulse with (a) Ne = 105 , (b) Ne = 106 and (c)
Ne = 107 electrons and (d) for an elliptical laser pulse with Ne = 106 electrons.

In conclusion, we have derived a mean field formalism to simulate the propagation of

86

Gaussian
Gaussian
Gaussian
Elliptical

-

Ne
Ne
Ne
Ne

= 105
= 106
= 107
= 106

zf ocus
(mm)
1088.54
1088.62
1088.57
1088.61

√

σT
(µm)
0.98
0.92
0.95
0.46

√

σz δt
(µm) (fs)
0.76 13.1
2.27 39.2
16.6 288
2.37 41.1

Table 3.2: Focal spot position and corresponding pulse width for varying number of electrons
Ne and transverse photoexciting laser profile.
an electron pulse through a microscope column with different optical elements. By comparing the time evolution of the pulse described with the Analytic Gaussian model to our
N-particle simulations, we have shown good agreement for extraction conditions under the
virtual cathode limit. Following a discussion of the effects of the different optical elements
composing the microscope column, we have built a model to describe our UEM system and
given quantitative predictions of the achievable spatial and temporal resolutions. Given the
limitations of the Analytic Gaussian model, these predictions should be viewed as best case
results, since aberrations and deviations of the pulse from the assumed self similar Gaussian
profile were neglected. Nevertheless, due to the simplicity of the model and the resulting
ease of computation, we believe these results can offer a valuable tool in understanding the
achievable resolution in UEM systems and allow one to search optimal configurations over a
vast range of parameters. Other simulation tools (such as the COSY-MLFMM code) could
then be used to analyze the most promising configurations in more detail.

87

Chapter 4
Microscopic origin of the properties of
TaS2
The final part of this work will focus on an example of a material, Tantalum Disulphide
(TaS2 ), that can be studied with the Ultrafast Electron Microscope described in the previous
chapters. To do this, we first give a brief introduction to correlated electron systems with an
emphasis on Charge Density Wave (CDW) insulators, followed by a survey of the existing
experimental and theoretical results on TaS2 . We then proceed to illustrate the simulation
methods used and conclude with a discussion of the results.

4.1

Physics of correlated electron systems

Materials with partially filled d- and f-orbitals have received considerable interest over the
last decades due to the challenges they present. The localized nature of these orbitals can
lead to the opening of a gap in the ground state in systems that would otherwise be metallic.
Moreover, varying the temperature (or other macroscopic parameters such as pressure or
88

doping), can lead to a transition from the insulating to the metallic state with a change in
resistivity of several orders of magnitude. Over the last 30 years explaining these interesting
features has been a major challenge to theory and experiment alike and has led to the
development of many novel techniques (for an introduction to the subject see [73] and [27]
and references therein).
The unusual behaviors seen in these so-called correlated electron systems represent a
many-body problem that arises from the interaction of electronic and structural degrees of
freedom. Based on the driving mechanism of gap opening, it is possible to define three
categories of insulators (Fig. 4.1). The first category is the so-called Mott insulator, in
which the electron-electron Coulomb repulsion drives a gap opening at the Fermi energy EF
and leads to localization of the valence electrons on the atomic sites. The second type of
insulator is the excitonic insulator, resulting from unscreened electron-hole interactions.
In these materials a Charge Density Wave (CDW) is formed together with a small distortion
of the lattice that is caused by a finite electron-phonon coupling. Due to an avoided crossing
of valence and conduction bands, a gap opens at EF that stabilizes this distortion. The
third type of insulating behavior is the result of electron-phonon interactions and is called
the Peierls insulator. In this case a gap opens as a result of the simultaneous formation
of a CDW and a periodic lattice distortion with a wavelength that is determined by the
periodicity of the Fermi surface. The energy of the occupied levels is lowered by an amount
that is greater than the Coulomb repulsion resulting from the distortion.
The term Charge Density Wave introduced here indicates a periodic modulation of the
ground state electron density. This phenomenon was first described by Peierls in 1955
[74] when he noted that a one dimensional metal is not stable at low temperature. A
periodic modulation of the charge density, accompanied by a periodic lattice distortion due
89

Figure 4.1: Three classes of insulators and their respective characteristic response times. (a)
Mott insulator (b) Excitonic insulator (c) Peierls insulator. In each case the gray shading
and dashed lines indicate the metallic state and the solid lines and colored shading show the
resulting insulating state both in real and momentum space. (d) Timescales corresponding
to probing these processes with a pump-probe experiment. [27]
to electron-phonon interactions, leads to an insulating ground state with an energy that is
lower than that of the metallic state. The size of the insulating gap depends on the amplitude
of the lattice distortion, while the filling of the levels in the system determines the periodicity.
For more details on CDW formation, we refer the reader to the classic text by Gruner [75].
Of course real materials present competing electron-electron and electron-phonon interactions and thus rarely fit in a single one of the categories described above. However, as shown
in the cartoon of Fig. 4.1(d), these interactions can be separated based on the temporal
response of the system to perturbation with a pump pulse. Electron hopping processes are
typically very fast and act on timescales of τ <10fs. The buildup of charge screening that
results in the formation of excitonic behavior is a slower process characterized by time constants τ =10-100fs. The excitation of phonon modes leading to charge oscillations requires a
transfer of energy from the excited electrons to the lattice and therefore takes place at later
times, with τ >100fs. Based on this separation of time scales, pump-probe measurements
have become a key technique to investigate these systems, as measuring the relaxation of

90

the system to equilibrium after perturbation with the pump pulse allows decoupling of these
degrees of freedom.

4.2

Structural and electronic properties of TaS2

For the present work we chose to investigate Tantalum Disulphide as this material is thought
to be a combined Mott-Peierls insulator exhibiting strong electronic localization, CDW formation and lattice distortion. TaS2 is a member of the transition metal dichalcogenides
group and presents a 2 dimensional layered structure, as seen in Fig. 4.2, with weak Van
der Waals interactions between the layers and covalent bonds between atoms in the same
layer [76,77]. The stacking of the layers is thought to be disordered [78–80]. The charge den√ √
sity wave distortion leads to the formation of a commensurate superlattice of size 13x 13
and rotated by 13.9◦ with respect to the original 1x1 unit cell for temperatures T <180 K [81].
In this “Star of David” arrangement, shown in Fig. 4.3, we can distinguish three types of
inequivalent atoms: one central atom which is not shifted with respect to the original lattice (site a), its six nearest (b sites) and next-nearest neighbors (c sites) which undergo a
contraction towards the central atom forming the clusters drawn in the figure. This in-plane
displacement is accompanied by an out-of-the-plane swelling of the sulfur layers [1].
The CDW distortion alone is not sufficient to fully open a gap in TaS2 , as only twelve of
the thirteen electrons of the reconstructed cell occupy levels below the distortion-generated
gap. Since Angle Resolved Photoemission Spectroscopy (ARPES) results [81–85] and resistivity data [83,86] show this material to be an insulator, some further interaction is needed to
explain these results. Conventionally it is thought that a gap opens as a consequence of the
splitting of the band corresponding to the thirteenth electron, localized on the central atom,

91

into upper and lower Hubbard bands due to Coulomb repulsion leading to the development
of a Mott insulating state [28]. This explanation for the origin of the gap is not entirely
convincing as some photoelectron spectroscopy (PES) investigations [87], tight-binding [81]
and DFT calculations [1, 88, 89] show that the system retains a finite density of states at the
Fermi energy. An alternative explanation suggests that the existence of a gap could be due
to Anderson localization caused by disorder [89]. The origin and nature of this insulating
gap is one of the topics that will be investigated in the present work.

Figure 4.2: TaS2 structure corresponding to the primitive unit cell (undistorted metallic
phase): (a) 3D view, (b) Top view and (c) Side view. Ta atoms are drawn in blue, S atoms
in red. The black lines delimit the primitive 1x1 unit cell.

Increasing the temperature leads to breaking of the commensurate CDW (C-CDW) and
the insulating state, and between T = 190K and T = 350K the system exhibits a nearly
commensurate CDW (NC-CDW) corresponding to a temperature dependent rotation of
the CDW vector from the commensurate value of 13.9◦ to about 11◦ and a simultaneous
formation of domains [77,90,91] with a size inversely dependent on temperature. For 350K <
T < 550K, we find an incommensurate CDW (IC-CDW) in which the CDW vector is aligned
with the underlying lattice and no domain structure is visible. For T > 550K the system
92

Figure 4.3: Positions of the Ta atoms corresponding to the reconstructed “Star of David”
arrangement in a 2x2 unit cell. The undistorted positions are shown in yellow, while the
reconstructed lattice and bonds resulting from the contraction associated with CDW formation are drawn in blue. The labels correspond to the three inequivalent atom types of the
distorted cell.
becomes metallic and reverts to a monoclinic undistorted lattice. Figure 4.4 shows resistivity
as a function of temperature and permits a clear identification of the phases discussed. Also
shown is the trigonal phase obtained upon heating between T = 200K and T = 280K [77].

Figure 4.4: Resistivity as a function of temperature in 1T-TaS2 . The phases shown are: CCDW below 190K, NC-CDW below 350K, IC-CDW below 550K and metallic for T>550K.
Also shown is the trigonal phase present during the heating cycle between 200-280K. [28]

93

Both pressure [28, 92, 93] and doping [94] can be used to modify the boundaries between
these phases and even induce the formation of a superconducting state. How this phase is
formed and whether it coexists with the CDW in real space or in reciprocal space is still an
unanswered question in the field and is outside of the scope of the present work.
As we have seen, experimental techniques that probe out of equilibrium conditions, such
as pump-probe measurements, can untangle the different timescales involved in the interactions between electronic and phononic degrees of freedom and offer an excellent tool to
better understand the complex phase diagram present in TaS2 . By performing time resolved
ARPES measurements, Perfetti et al. [95, 96] showed the ultrafast collapse of the electronic
gap after photoexcitation due to the creation of a hot electron distribution, and its successive revival after 680 fs. The nuclear and phononic degrees of freedom, acting on longer time
scales due to the coupling to the electronic ones, lead to periodic oscillations of the system
at long times (t>1-2ps) due to excited phonon modes. A qualitative explanation of these
results using a single band Hubbard model in the context of DFT+DMFT was advanced by
Freericks et al. [97], but several discrepancies such as a strongly insulating behavior absent
in the experiments, limit the applicability of these results.
More recently, time resolved x-ray photoemission spectroscopy (TR-XPS) has been used
to probe the melting of the CDW at specific atomic sites [27, 86] and the transition from CCDW to NC-CDW. The authors find an initial ultrafast subpicosecond reduction of the CDW
peak and its successive partial recovery with a time constant of ∼900 fs. This seems to suggest
that the system relaxes into a state with an inhomogeneous distribution of undistorted
(metallic) islands inserted into a CDW background, instead of recovering the superstructure
network typical of the equilibrium NC-CDW.
Recent experiments were also conducted by Prof. Ruan’s group at Michigan State (cur94

rently unpublished) to understand the laser fluence and wavelength dependence of the transition from NC-CDW to IC-CDW in TaS2 . By measuring both the CDW peak amplitude
and rotation angle, they observe a transition dependent on the incident photon density,
rather than on energy density as would be expected. This suggests that the transition is not
thermodynamically driven but rather driven by a direct distortion of the CDW.

4.3

Methods

To better understand the ground and excited state properties of TaS2 , we performed ab-initio
Density Functional Theory (DFT) calculations of the electronic and phononic structure both
in the commensurate CDW phase and in the undistorted metallic phase using the Vienna Abinitio Simulation Package (VASP) [98–101]. A description of the basic equations governing
the DFT theory is given in appendix 5. For the calculations presented in the rest of this
chapter, the GGA-PBE [102] exchange correlation functional was used and plane waves with
the projector augmented wave method [103] were used as the basis set.
The lattice constants and atom positions for the undistorted metallic phase (1 Ta atom
√ √
and 2 S atoms in the unit cell) were taken from Spijkerman et al. [90]. The 13x 13 CCDW supercell (13 Ta atoms and 26 S atoms) with a 13.9◦ rotation was then built from
this primitive unit cell using the bond contractions reported by Smith et al. [81]. For each
structure, the first step was to determine the appropriate energy cutoff and k-point mesh
such that the total (ground state) energy was converged to 1meV/atom. The results of these
checks are show in Fig. 4.5 for the undistorted cell and in Fig. 4.6 for the C-CDW supercell.
Based on this data, we choose to use an 11x11x6 mesh with an energy cutoff Ecut = 300eV
for the metallic system and, for the C-CDW system, a 4x4x8 mesh with Ecut = 350eV.

95

Ground state energy (eV)

-23.14
-23.145
-23.15
-23.155
-23.16
-23.165
-23.17
-23.175
-23.18
200 250 300 350 400 450 500
Energy cutoff (eV)

5

10
15
k-point mesh

20

Figure 4.5: Convergence of the ground state energy for the metallic phase as a function of
energy cutoff and k point mesh size. The dashed lines correspond to the acceptable range of
values that gives a precision of 1meV/atom. Note that due to the size of the reciprocal unit
cell, the Nth k-point mesh corresponds to a mesh of size N in the kx and ky directions and
N/2 in the kz direction.

Ground state energy (eV)

-301
Ecutoff = 350 eV
Ecutoff = 300 eV

-301.1
-301.2
-301.3
-301.4
-301.5
-301.6
300

350 400 450 500
Energy cutoff (eV)

550 1

2

3

4
5
6
k-point mesh

7

8

9

Figure 4.6: Convergence of the ground state energy for the distorted phase as a function of
energy cutoff and k point mesh size for varying energy cutoff. The dashed lines correspond
to the acceptable range of values that gives a precision of 1meV/atom. Note that due to the
size of the reciprocal unit cell, the Nth k-point mesh corresponds to a mesh of size N in the
kx and ky directions and 2 · N in the kz direction.
Next, the positions of the atoms in the C-CDW structure were relaxed by performing
self consistent DFT runs with a force calculation. The conjugate gradient algorithm was
first used to relax the structure until the force on each atom was less than 0.1eV/˚
A and
the total and band structure energy changes between steps in the self-consistent electronic
96

loop were less than 10−2 eV. This was followed by a more stringent run in which a residual
minimization scheme - direct inversion in the iterative subspace (RMM-DIIS) algorithm [104]
was used to relax the forces to less than 0.01eV/˚
A, with total and band structure energy
changes ∆E < 10−5 eV. The size of the unit cell in all three dimensions was allowed to
vary in addition to the atomic positions. To get accurate values for the stress tensor used to
calculate the forces on the atoms, it is necessary to increase the energy cutoff when performing
structural relaxations, so a cutoff equal to 1.3 · Ecut was used, following the recommendation
in the VASP documentation. The optimized unit cell with the respective atomic positions
was then used in all the following calculations.
To simulate the metallic state originating from pumping with an ultrafast laser, the
optimized unit cell for this phase was built using the lattice constants obtained from relaxing
the C-CDW structure, with the atoms restored to the initial high symmetry positions. Using
the symmetry present in the system, the unit cell of the metallic phase can be reduced to a
cell with 1 Ta atom and 2 S atoms.
For both the metallic and the C-CDW phase, we calculated the band structure and the
density of states (DOS); for the latter a denser k-point mesh was used equal to, respectively,
21x21x11 (metallic phase) and 8x8x16 (C-CDW phase). The effects of the spin-orbit coupling were also investigated, as previous results from a tight-binding calculation [88] suggest
that they play an important role in TaS2 . To see whether a Mott gap would open in the
C-CDW phase, DFT+U [105] was used with varying values of on-site Hubbard interaction,
U = 0 − 5eV. Phonon spectra for both the metallic and C-CDW phases were also calculated
using density functional perturbation theory [106] in VASP and the results were post processed using the program Phonopy; for the metallic phase a 3x3x3 unit cell was used in this
calculation with phonon wave vectors sampled on a 12x12x6 mesh. For the C-CDW phase
97

we used a 1x1x1 unit cell and an 8x8x16 mesh. Finally, to investigate the optical absorption
of these materials, we computed the real and imaginary parts of the dielectric constants in
the independent particle approximation [107].

4.4

Ground state of TaS2

We begin by presenting our efforts in understanding the ground state of TaS2 and the origin
of the gap that is seen in the experimental data. To do so we calculate the band structure and
density of states and investigate whether spin orbit coupling or Hubbard splitting could lead
to an insulating state in this material. Following this, we will discuss the role of phonons in
thermalizing the hot electron distribution generated by the pump laser used in experiments
and conclude with a brief discussion of the optical properties of TaS2 .

4.4.1

Band structure and DOS

After performing the relaxation of the structural parameters (unit cell vectors and atomic
positions) for the unit cell corresponding to the C-CDW phase, as described in the previous
section, we find the atomic ground state configuration shown in Fig. 4.7, with lattice parameters indicated in Table 4.1. As expected, we observe a buckling of the S layers both above
and below the Ta atomic plane with a variation in height of the S atoms equal to ∆z = 0.2 ˚
A.
The lattice parameters obtained show good agreement with experimental data, with a and b
showing a 0.5 % difference from the experimental value. The vertical c axis is elongated by
about 13% compared to the experimental result due to the weak Van der Waals interactions
between the layers that is underestimated by DFT. The first Brillouin zone corresponding
to the C-CDW tetragonal unit cell is shown in Fig. 4.8, with the high symmetry points that

98

will be discussed in the following text labeled.

Figure 4.7: Top (a) and side (b) views of the relaxed C-CDW structure, with the black line
indicating the boundary of the unit cell. Ta atoms are drawn in blue, while S atoms are in
red.
Theory
Exp
a (˚
A) 12.192 12.129
b (˚
A) 12.192 12.129
˚
c (A)
6.677 5.889
Table 4.1: Theoretical and experimental (from [1]) values for the parameters of the C-CDW
unit cell

Figure 4.8: First Brillouin zone corresponding to the C-CDW TaS2 unit cell. The labels
represent the conventional high symmetry points. [29]
Our calculated band structure and density of states (DOS) are reported in Fig. 4.9
and agree with similar data previously published in the literature [1]. The uppermost half
99

occupied band shows a weak dispersion in plane (Γ-M-K-Γ), being localized at energies of
about −0.1 eV and a strong dispersion in the out of plane direction (A-Γ, Γ-L, H-Γ) that
results in a crossing of the Fermi energy at positions close to the Γ point. This band is
mostly due to the Ta dz2 orbitals. Below this, at energies −0.2 eV to −1.0 eV, we find the
six sub-bands corresponding to the remaining d orbitals. These bands are somewhat close
at the Γ point (−0.3 eV to −0.5 eV) but elsewhere display a large spread in energies. The
DOS plot shows the presence of states at the Fermi energy, mostly due to the dz2 band on
the central Ta atom, so that the system is not fully gapped.

√13 x 

√13 x 1 - C-CDW
1
0.5
E (eV)

0
-0.5
-1

total
Ta dz2
Ta d

-1.5
-2
A

Γ

MK

Γ

L H

Γ

dos (a.u.)

Figure 4.9: Band structure along high symmetry lines and density of states for the C-CDW
phase. The contribution to the DOS due to the Ta d orbitals and that due to only the dz2
is also plotted.

In contrast the band structure displayed by the metallic case, Fig. 4.10, has an uppermost
band which is localized in the vertical A-Γ direction and strongly dispersive along the other
directions considered. From this we can deduce that the CDW distortion acts in localizing
the dz2 orbital corresponding to the central Ta atom in plane with the swelling of the S
planes caused by the consequent delocalization in the vertical direction. The DOS plot for
the metallic phase shows that the density at the Fermi level includes a significant contribution
100

due to the d orbitals which, unlike what seen in the C-CDW phase, is not fully originated
by the dz2 orbital.

1 x 1 x 1 - Metallic
3
2

total
Ta dz2
Ta d

E (eV)

1
0
-1
-2
-3
-4
A Γ

M K

Γ

L H

Γ

dos (a.u.)

Figure 4.10: Band structure along high symmetry lines and density of states for the metallic
phase. The contribution to the DOS due to the Ta d orbitals and that due to only the dz2
is also plotted.

Comparing the ground state energies (Table 4.2), we see that the distorted C-CDW is
indeed more stable at low temperatures compared to the metallic phase and that, even
though it is not completely gapped, the density of states at the Fermi energy is reduced
fivefold.
E0 (eV/TaS2 )
n(EF ) (states/eV/TaS2 )

C-CDW
-23.190
0.290

Metallic
-23.145
1.475

Table 4.2: Ground state energy E0 and density of states at the Fermi energy n(EF ) in the
C-CDW and metallic phases.

4.4.2

Role of spin-orbit coupling

Having investigated the ground states of the metallic and C-CDW phases, we now turn our
attention to the effect of the spin orbit coupling on the energy levels of TaS2 . A previous 2D
101

tight binding calculation [88] suggested that this coupling could split the band corresponding
to the dxy , dx2 −y2 orbitals on the centralized atom, generating a narrow band of states at
the Fermi level. DFT calculations on TaS2 were performed in the presence of spin-orbit
coupling [97] but its role on the band structure and density of states was not fully explored.
Figure 4.11 presents our results for the C-CDW phase both with and without spin orbit
coupling. Bands below EF (doubly occupied) and those above EF (unoccupied) are slightly
perturbed by introducing this coupling, while the upper half filled band is relatively insensitive to it. The density of states also shows a slight perturbation as a result but no significant
changes are observed. To better understand the effect of the spin orbit coupling, we focus
on the density of states around EF in Fig. 4.12. The atom projected densities [panel (a)]
confirm our previous observation that the DOS at the Fermi level is mainly localized on the
central Ta atom (a-site - see Fig. 4.3). This remains the case even with the introduction
of the SO-coupling, and its main effect is seen in a shift the first peak in the density above
EF from 0.2 eV to 0.25 eV. The orbital character of the density of states [panel (b)] shows
that at EF the states are still mostly of dz2 character. We conclude that the splitting of the
antibonding dx2 −y2 and dxy bands seen in the 2D tight binding calculations is suppressed in
3 dimensions, where the uppermost band retains its dz2 character with a strong out of plane
dispersion.
We also report the band structure and DOS for the metallic case in the presence of SOcoupling in Fig. 4.13. In this case, in agreement with the previous calculation by Rossnagel
et al. [88], the SO-coupling is seen to lift the degeneracy of the dz2 and dxy , dx2 −y2 at the Γ
point and induce an energy splitting of 0.3 eV.

102

C-CDW
C-CDW (+ SO coupling)

1

E (eV)

0.5
0

-0.5
-1
A

MK

Γ

Γ

L H

dos (a.u.)

Γ

Figure 4.11: Band structure along high symmetry lines and density of states for the C-CDW
phase both with (black) and without (red) spin orbit coupling.

(a) Atom type
dos per atom (a.u.)

a
b
c

-0.4

(b) Orbital type
dz2
dx2-y2+dxy
dxz+dyz

DFT
DFT+SO

-0.2

0

0.2

0.4

E (eV)

-0.4

-0.2

0

0.2

0.4

E (eV)

Figure 4.12: Detail of the density of states for the C-CDW phase both with (continuous
lines) and without (dashed lines) spin orbit coupling. (a) Atom projected densities. (b)
Orbital projected densities.

4.4.3

Role of Hubbard repulsion

Since spin orbit coupling acts as a perturbation and does not play a significant role in
determining the ground state properties of TaS2 , we now turn our attention to understanding
the effect that an on-site Hubbard repulsion has on the ground state of the C-CDW phase

103

Metallic
Metallic (+ SO coupling)

3
2
E (eV)

1
0
-1
-2
-3
A Γ

M K

Γ

L H

Γ

dos (a.u.)

Figure 4.13: Band structure along high symmetry lines and density of states for the metallic
phase both with (black) and without (red) spin orbit coupling.
and whether there is evidence of a Mott gap opening. The results for varying values of the
electron-electron interaction U are plotted in Fig. 4.14. The first thing we observe is that
there is no evidence of a gap opening even at the highest value of U considered (U=5 eV).
While this is in agreement with the data presented in a preprint by Darancet et al. [89],
previous DFT+DMFT calculations [96, 108, 109] showed a gap opening by using a single
band Hubbard model. The origin of the discrepancy it most likely due to the simplified
model used in these papers that neglects the effect of the neighboring occupied d orbitals
and of the coupling to any unoccupied states.

4.4.4

Role of stacking disorder

Since we do not see evidence of the formation of a Mott gap due to Hubbard repulsion,
our results suggests that TaS2 is an in-plane insulator but presents metallic states in the
out of plane direction and that the experimentally seen gap is originated in the disordered
stacking sequence present in this material. We now proceed to test this hypothesis by using

104

0.4
dos (a.u.)

E (eV)

0.2
0
-0.2

U=0.0 eV
U=1.0 eV
U=2.3 eV
U=5.0 eV

-0.4
A

Γ MK Γ

LH

Γ

-0.4 -0.2

0
0.2
E (eV)

0.4

Figure 4.14: Band structure along high symmetry lines and density of states for the C-CDW
phase in the DFT+U approximation with varying values of the on-site electron-electron
interaction U
the stacking sequence for TaS2 proposed in [80], in which the authors find good agreement
with the structure factor calculated from the X-ray data in [78, 79] and suggest that the
stacking is of type ACi ACj ..., where A indicates stacking of layers with the central a-type
atoms lying on top of each other, while Ci indicates configurations where the central a-type
atom lies on top of a c-type atom (see Fig. 4.3 for the definition of a and c sites). There
are six possible Ci type stackings due to the six c atoms present in each unit cell and in the
discussion that follows we will collectively indicate them simply as C stacking. The order in
which these C stacking units repeat is random.
As density functional theory calculations scale approximately as O(N 3 ), where N is the
number of atoms in the system, we are limited in how many layers we can consider. Therefore
we limit ourselves to the three cases: systems with A stacking (which is what we have used
in the calculations presented so far and corresponds to a configuration with no disorder), C
stacking and AC stacking; the unit cells used for each system are shown in Fig. 4.15. Note
that we use periodic boundary conditions in all three directions, which is the reason for using
4 layers in the unit cell for the AC stacking configuration.
105

Figure 4.15: Units cells for the stacking sequences considered: (a) A stacking, (b) C stacking
and (c) AC stacking.
For each stacking configuration we calculate the band structure and density of states and
present our results in Fig. 4.16. In the absence of stacking faults, as we already showed in
the previous sections, the system presents a finite density of states at the Fermi energy and
no evidence of gap opening. Considering two layers with C stacking, we see that each central
a-atom contributes one weakly dispersing band around the Fermi energy. While the broad
feature seen in the A stacking for an energy range of −0.25 eV to 0.125 eV is narrowed, the
two bands in the C stacking are not fully separated in energy so that no gap opens in this case
either. For the AC stacking there are four bands in the energy range −0.25 eV to 0.25 eV,
and due to their mutual repulsion, we observe a gap in the density of states.
106

(a)

0.5

E (eV)

0.25
0
-0.25
-0.5
A

Γ

MK

Γ

L H

Γ

A stacking
(b)

dos (a.u.)
C stacking

0.5

E (eV)

0.25
0

-0.25
-0.5
A

Γ

M K

Γ

L H

Γ

A stacking
(c)

dos (a.u.)
AC stacking

0.5

E (eV)

0.25
0

-0.25
-0.5
A Γ

M K

Γ

L H

Γ

dos (a.u.)

Figure 4.16: Band structure along high symmetry lines and density of states for the C-CDW
phase with different stacking: (a) A stacking, (b) C stacking and (c) AC stacking. The
last two panels also show a comparison to the DOS obtained with the A stacking (stacking
without disorder).
To check whether the AC stacking is indeed a possible configuration for the system, we
consider the ground state energy for each configuration and find values for the A and AC
stackings of −301.465 eV/layer and −301.455 eV/layer, respectively, while the C stacking
107

presents a ground state energy of −301.402 eV/layer. This means that for T >118K the A
and AC stacking sequences differ by less than the thermal energy. For the C stacking to be
within kT of the other two, a temperature of T =730K is required. These results confirm
that the ground state of TaS2 does indeed exhibit random stacking and show that contrary
to the conventional wisdom, the presence of a gap in the ground state of this material might
be due to stacking disorder rather than Mott physics.

4.5

Phonons

Having investigated the ground state of TaS2 both in the C-CDW and metallic phases, we
now want to understand how the system relaxes after being perturbed with a pump laser
pulse. tr-ARPES data [95] of the C-CDW phase shows a breakdown and successive recovery
of the metal insulator gap 680 fs after photo excitation and a subsequent long lived oscillation
of the energy of the system, with a frequency of about 2.5 THz (see Fig. 4 in [95]). We want
to explain how the excess energy of the hot electrons is transfered to the lattice and what
determines the temporal evolution seen in the experiments.
By using density functional perturbation theory [106], we calculate the phonon band
structure of both the C-CDW and the metallic phases and present it in Fig. 4.17. The
acoustic mode dispersion curves for the metallic state [panel (a)] show characteristic instabilities expected from a system with CDW formation [75], and the wave vector associated
with this instability coincides with the CDW wave vector. The negative values correspond
to imaginary frequencies representing unstable modes. The optical branches located above
20 meV are mostly due to vibrations of the lighter S atoms. Overall these results show a
mostly in-plane dispersion and localized phonon bands in the vertical direction (Γ-A) and

108

agree with a previous calculation on phonons in the metallic state of TaS2 reported in the
literature [92].

(b) C-CDW

50

10

40

8

30

E (meV)

E (meV)

(a) Metallic

20

6
4

10

2

0

0
A Γ

M

K

Γ

(c)
2

4

dos-per-TaS2 (a.u.)

0

A

frequency (THz)
6
8

M

Γ

10

K

Γ
12

C-CDW
Metallic
C-CDW - Ta
Metal - Ta

0

10

20

30
Energy (meV)

40

50

Figure 4.17: Calculated phonon band structure for the metallic (a) and C-CDW (b) phases.
Imaginary frequencies, corresponding to unstable modes are represented by negative values.
(c) Phonon density of states for the metallic (green) and C-CDW (black) phases. The
continuous line represents the total DOS, while the dashed line shows the contribution due
to the Tantalum atom vibrations.

On the other hand the phonon band structure for the C-CDW phase [panel (b)] displays
a wealth of phonon branches due to the numerous atoms present in the unit cell. Only the
lower ones are shown in Fig. 4.17(b); an upper band of modes in the range 20 meV to 50 meV
is present but it is not of interest for the current discussion as it is mostly due to motions

109

of the S atoms. Again we observe localized phonon bands in the vertical direction (Γ-A),
indicating soft phonon modes. There is also evidence of avoided crossing in the in plane
(Γ-M-K-Γ) direction, which points to strong interactions between phonon modes leading to
mode mixing. In particular in the Γ-M segment, at about 2/3 of the way to M, we see an
avoided crossing of the LA and LO bands, which presents itself again about halfway between
M-K and at K. Other avoided crossings are also present in the higher energy optical modes.
The phonon density of states [Fig. 4.17(c)] shows weak mixing between the Ta and S
atom vibrations, with a clear separation of the frequency ranges due to each one. Since
the distortion of the CDW is due to motions of the Ta atoms, we focus on the low energy
range and present in Fig. 4.18, in addition to the total phonon DOS, the projection of
the DOS onto the atomic sites in the reconstructed unit cell. In the C-CDW phase, we
observe numerous phonon peaks between 8 meV to 16 meV, with the principal one at 13 meV
(corresponding to a frequency of 3.1 THz) due to a vibration of the central atom (a-site). In
the metallic phase, on the other hand, the main mode is localized at 10 meV (2.3 THz), with
a lower energy broad secondary peak at 6 meV (1.5 THz).
As the system is in the metallic state following the pump laser excitation, the phonon
DOS for this state will determine what modes are most likely to be excited to dissipate
the excess energy. Our data shows the mode at 2.3 THz to be the most likely candidate,
in striking agreement with the experimentally measure oscillation frequencies reported by
Perfetti (see Fig. 4 in [95]). In the C-CDW phase, this vibrational frequency corresponds to
a superposition of modes involving all Ta atoms in the reconstructed cell. These modes then
can dissipate their energy through interaction with lower energy modes, seen in the avoided
crossings discussed in relation to the band diagram (Fig. 4.17).

110

0

frequency (THz)
2

1

3

4

dos-per-TaS2 (a.u.)

C-CDW
C-CDW - Ta A
C-CDW - Ta B
C-CDW - Ta C
Metallic

0

2

4

6

8
10
12
Energy (meV)

14

16

18

Figure 4.18: Detail of the phonon density of states per TaS2 for the C-CDW phase in the
range 0 meV to 18 meV (corresponding to 0 THz to 4 THz). The red lines indicate the site
projected DOS.

4.6

Optical absorption

We conclude the discussion by presenting the results relative to the optical absorption profile
of this system, with the goal of understanding the optimal conditions to excite this material
in the experimental setup. The real and imaginary parts of the dielectric function are
computed in the independent particle approximation [107] in which classical depolarization
effects (local-field effects), electron-hole and electron-electron interactions are not included.
The real part

1 (ω)

of the dielectric function describes the polarization induced by the field

and is therefore related to electron dispersion inside the material, while the complex part
2 (ω)

describes the absorption properties. Other optical functions that can be calculated

from the dielectric function are:

111

• the complex refractive index n(ω) + ik(ω):

1

n=

= n2 − k 2 ,
2
1

+ 22 +
2

1

,

2

= 2nk,

(4.1)
2
1

k=

+ 22 −
2

1

,

(4.2)

• the optical absorption coefficient α(ω):

α=

4π
ωk(ω),
hc

(4.3)

• the electron energy-loss function:

−

1
(ω)

=

2
2
1

+

2
2

,

(4.4)

• the reflectance R(ω):
R=

(n − 1)2 + k 2
,
(n + 1)2 + k 2

(4.5)

Note that whenever necessary, we have dropped the frequency dependence from the previous
equations for simplicity.
Our results for the in-plane direction are summarized in Fig. 4.19, where the first thing
we note is a remarkable agreement between our simulated data and the experimental measurements found in the literature [see Fig. 10(b)-(d) of [76]]. The energies of the two main
peaks of the imaginary part of the dielectric function in the C-CDW phase [panel (b)], corresponding to intraband transitions, are well described by our simulations and match with
the experimental results. The shoulder present in the peak at 1.5 eV is also captured in our

112

calculations. The main difference present is an overestimation of the value of the imaginary
part and an underestimation the value of the real part. This has the effect of introducing
additional zero-crossings of

1

at energies < 5eV that are not seen in the experimental data

and is most likely due to the absence of impurities or defects in our simulated system. The
high positive value of

1 (0)

and the absence of a metallic behavior [ 1 (0) < 0], suggests that

TaS2 has a small optical band gap.
The absorption coefficient α [panel (c)], is also accurately described, as the shapes of
the secondary peak at 2.5 eV and of the broad feature between 5 eV and 20 eV are captured.
Similar considerations hold for the energy-loss function in panel (d), where the small peak
at 2.5 eV due to interband transitions is clearly captured, in addition to the plasmonic peak
seen at 20 eV [76]. The strength of this peak is increased by an order of magnitude compared
with the experimental data, which again can be explained by the absence of imperfections
in the simulated material.
The results for the metallic phase are also reported in Fig. 4.19, but no experimental
verification has been found in the literature for this data. Due to the smaller size of the unit
cell used, the plots are overall noisier. The absorption coefficient and energy-loss function
retain a behavior similar to that seen in the C-CDW phase, even though the details of the
dielectric function are modified with the metallic phase displaying a peak in

2

at an energy

of 0.75 eV, absent in the C-CDW state.
Given the good agreement with experiment, we conclude that the optical transitions in
TaS2 are mostly single particle excitations from occupied valence states to unoccupied conduction states, with electron-electron and electron-hole interactions and polarization effects
playing a minor role. Some differences were seen in the magnitude of the optical properties
studied, which are due to imperfections in the sample used in the experiments that are not
113

(a) Metallic
35
30
25
20
15
10
5
0
-5

(b) C-CDW
50

30
20
10
0

45
30
15
0

40
30
0

1

2

20

3
ε1
ε2

0

10

2

3
ε1
ε2

0
0

5

10 15 20
Energy (eV)

25

0

5

(c)

10 15 20
Energy (eV)

25

(d)

2.5

3
Metallic
C-CDW

2

Metallic
C-CDW

2.5
Im(-1/ε)

α (x102 1/µm-1)

1

1.5
1
0.5

2
1.5
1
0.5

0

0
0

5

10 15 20
Energy (eV)

25

0

5

10 15 20
Energy (eV)

25

Figure 4.19: In-plane optical parameters for TaS2 : real ( 1 ) and imaginary ( 2 ) parts of the
dielectric constant for the metallic (a) and C-CDW (b) phases. The insets show a detail
of the dielectric constant in the energy range relevant for pump-probe experiments. (c)
Absorption coefficient (d) Energy-loss function
captured in the simulations. In general we expect the inclusion of local fields to not play
a significant role in the in-plane components as the system is mostly homogeneous, while
it could have an effect if one is interested in the out of plane optical properties [110]. The
peak seen in

2

at 1.0 eV suggests that exciting TaS2 with a wavelength corresponding to this

energy, would give the highest absorption and allow probing the material with a lower laser
power and higher efficiency.

114

For completeness we also report the predicted out of plane optical parameters for TaS2
in Fig. 4.20. We note the overall lower value of the imaginary part of the dielectric constant, indicating, as expected, that the material is more absorbing in plane. The absorption
coefficient and energy loss function show features very similar to those already discussed for
the in-plane direction.

(a) Metallic

(b) C-CDW

15

15
ε1
ε2

10

ε1
ε2

10

5

5

0

0

-5

-5
0

5

10 15 20
Energy (eV)

25

0

5

(c)

25

(d)

2.5

3
Metallic
C-CDW

2

Metallic
C-CDW

2.5
Im(-1/ε)

α (x102 1/µm-1)

10 15 20
Energy (eV)

1.5
1
0.5

2
1.5
1
0.5

0

0
0

5

10 15 20
Energy (eV)

25

0

5

10 15 20
Energy (eV)

25

Figure 4.20: Out of plane optical parameters for TaS2 : real ( 1 ) and imaginary ( 2 ) parts
of the dielectric constant for the metallic (a) and C-CDW (b) phases. The insets show a
detail of the dielectric constant in the energy range relevant for pump-probe experiments.
(c) Absorption coefficient (d) Energy-loss function

115

Chapter 5
Conclusions
The results presented in this thesis show how one can approach Ultrafast Science both from
the perspective of instrumentation development, as the achievable spatial and temporal
resolutions are fundamental quantities that determine what phenomena one can hope to
observe, and from the perspective of materials applications that can be studied with the
UEM system.
We started in chapter 2 by simulating the electron pulse photoemission process. By looking at the time dependent evolution of a single pulse and comparing it with experimental
results obtained using the shadow imaging technique, we confirmed the validity of our approach and saw that overall, for short timescales, the spatial profile of the pulse is dependent
on the transverse laser profile used to generate it, with velocities showing a turbulent flow
related to the initial thermal distribution and independent of extraction field. At later times,
as the pulse moves away from the surface, the pulse properties display a strong dependence
on the extraction field: for high fields the pulse becomes fully laminar, evolving into the
typical pancake shape, while for the lower extraction fields it retains some turbulent flow
and does not fully detach from the surface. By fitting the resulting transverse profiles with
116

Elliptical and Gaussian functions, we observed the surprising emergence of a 2D ellipsoidal
profile from the complex interplay of image- and space charge effects.
Our investigation of the optimal conditions for pulse generation in terms of laser fluence
and extraction field Fa for different experimental realizations showed evidence of virtual
cathode (VC) formation as the initial number of electrons was increased due the space charge
effect. By approximating the electron pulse as a sheet of charge, we derived a relationship
connecting the extraction field to the critical number of electrons required for the onset of
the VC.
We also analyzed the role of different photocathode geometries by observing the effect
that pinning of the image charge on the surface had on the onset of the VC limit and on the
final pulse parameters. In this case the positive image charge field provided a focusing effect,
increasing the fraction of electrons at the center of the pulse and preventing the pulse from
becoming fully pancake-like. Due to this both the transverse and longitudinal emittances
were seen to increase significantly once the virtual cathode limit was reached compared to
the case of radially expanding charge.
We showed that changing the shape of the transverse laser profile used to emit the
electrons also has a significant impact on the normalized emittance. Using a laser pulse with
a transverse elliptical profile gives the optimal pulse properties below the VC limit, while
once this limit is reached both the elliptical and uniform (top-hat) profiles present a sharp
scaling with the number of electrons, whereas the Gaussian profile retains a slightly less
unfavorable scaling. An increase in the extraction field can shift the onset of the VC regime
to higher numbers of electrons, thus increasing the favorable regime for the elliptical pulse.
The longitudinal emittance was observed to scale linearly with the number of electrons for
all laser profiles considered, with a weak dependence on the extraction field.
117

We conclude that optimal extraction conditions depend on a delicate balance between
the image charge, which disrupts the favorable internal self-fields of both the elliptical and
uniform cases, the number of electrons emitted and the extraction field, which controls the
onset of the virtual cathode limit. Below this limit the image charge has a weak effect with the
pulse properties primarily controlled by the laser parameters (transverse profile and fluence),
while above the VC limit the image charge strongly perturbs the expansion of the pulse with
a nonuniform dependence on the pulse profile. By calculating the brightness, coherence
length and energy spread for the pulses considered, we provided quantitative guidelines for
pulse generation.
In chapter 3 we built on our knowledge of the initial phases of pulse generation to simulate
the UEM system. We derived a mean field formalism to simulate the propagation of the
electron pulse through the microscope column, including the effects of the photoemission
gun, magnetic lenses, RF cavity and relativistic motion of the electrons. By comparing the
time evolution of the pulse described with this Analytic Gaussian model to our previous Nparticle simulations, we showed good agreement for extraction conditions under the virtual
cathode limit for both the Elliptical and Gaussian laser profiles. After analyzing the effects
of the different optical elements composing the microscope column, we built a model to
describe our whole UEM column and gave quantitative predictions of the achievable spatial
and temporal resolutions.
Moving forward, the next steps in the development of UEM systems will need to move
on two fronts. On the photoemission side, a more detailed model for the positive image
charge is needed and, since increasing the extraction field was seen to offer many advantages,
strategies to do this should be considered, including the addition of an RF cavity to the
photoemission gun, which would allow higher peak values of the electric field at the price
118

of added complexity in the system. On the other hand, to better understand the operation
of the whole UEM column, more accurate models and more sophisticated simulation tools
(such as the COSY-MLFMM N-particle code) are necessary to analyze in more detail the
configurations obtained with our simple Analytic Gaussian model. Given the limitations
of this model, the predictions we calculated should be viewed as best case results, since
aberrations and deviations of the pulse from the assumed self similar Gaussian profile were
neglected. A careful process of modeling the details of the lenses and cavity used in the
actual microscope is needed, so that details of their nonlinear effects can also be included in
the simulations.
On the applications side, in the last part of this thesis (chapter 4), we showed our work
in understanding the properties of TaS2 , a very promising material that can be studied with
the UEM microscope. With our density functional theory calculations of the electronic band
structure of 1T-TaS2 , we shone light on the ground state of this material and quantified
the effect of spin orbit coupling and Hubbard repulsion both in the metallic and in the
commensurate charge density wave (C-CDW) phases. Since our results showed that neither
of these interactions is sufficient to reproduce the insulating gap seen in experiment, we also
analyzed the effect of different stacking configurations of the TaS2 layers and found evidence
of gap opening for bilayers in the presence of disordered AC stacking. This points to a
deviation from the conventional picture of TaS2 as a Mott insulator, as the measured gap
could be due to the repulsion of the levels originating from interactions in the vertical stacking
direction. To better understand the processes involved in pump-probe measurements of this
material, we also looked at the phonon band structure and absorption spectrum and found
a remarkable agreement with the experimental measurements, which allowed us to point
out possible mechanisms for dissipation of the pump energy and recovery of the insulating
119

ground state.
Our understanding of stacking disorder in TaS2 opens the door to the investigation of
a variety of other layered materials (such as TaSe2 or TiS2 ) whose physics might also be
influenced by interlayer coupling, a mechanism that hasn’t been considered in the literature
so far. In addition to this, TaS2 itself still presents numerous unanswered questions such as
the emergence of superconductivity under pressure or doping, which might be better understood on the basis of our results. On a more general level, the role of disorder as a competing
mechanism in addition to Mott and Peierls physics in correlated electron materials is an open
question. Simulating disordered systems often requires significant computational resources,
but the inclusion of this effect might be necessary to explain some of the experimentally
observed material properties.

120

APPENDICES

121

Appendix A:
Generation of electron coordinates
This appendix gives details on how the initial x and y coordinates are generated in the
photoemission simulations for the three laser pulse shapes considered: Gaussian, Uniform
and Elliptical.

A.1

Generation of random numbers with an arbitrary distribution

The standard random number generator available in any programming language will provide
numbers with a uniform distribution in the interval (0, 1), but for many problems we wish to
use numbers generated at random from an arbitrary distribution. This section will outline
the method used to do so.
Starting with the given arbitrary normalized probability distribution, f (x), the first step
is to calculate the cumulative probability:
x

g(x) =

dx f (x ),

(A.1)

−∞

where the function g(x) represents the sum of probabilities from −∞ to x and has values
between (0,1), while x has values (−∞, ∞). We can now think of g(x) as a one-to-one
122

mapping from the interval of values of x (−∞, ∞) to (0,1). For this reason if we invert g(x)
by setting g(x) = y and calculating x = h(y), we have a function that maps the y values
belonging to the interval (0,1) into x. If we extract random numbers (y1 , y2 , ...yn ) from a
uniform distribution in the interval (0, 1), then the numbers (h(y1 ), h(y2 ), ...h(yn )) will be
random numbers from the original distribution f (x).
The tricky step in this procedure is inverting g(x), since for an arbitrary distribution it
might not be possible to find an analytic expression for its inverse. In that case numerical
methods would have to be used as an approximation.

A.2

Gaussian laser profile

In the case of a Gaussian profile, the distribution of coordinates is given by
x2 + y 2
1
exp −
f (x, y) =
2πσr2
2σr2

.

(A.2)

We begin by rewriting this distribution in terms of polar coordinates r, θ, where x = r cos θ,
y = r sin θ:
f (r, θ) =

r2
1
exp
−
2πσr2
2σr2

.

(A.3)

Since this distribution only depends on r, it is clear that we can simply choose θ from a
uniform distribution (0, 2π), while for r we will use the procedure presented in the previous
section.

123

We begin by calculating the cumulative probability distribution
2π

g(r ) =

r

dθ
r

2π

dθ

=

0

0
r

drr

=

drrf (r, θ)

(A.4)

r2
1
drr
exp − 2
2πσr2
2σr

(A.5)

0

0

0

1
r2
exp
−
σr2
2σr2

= 1 − exp −

r2
2σr2

(A.6)

.

(A.7)

Now we set g(r ) = α and solve this for r :

α = 1 − exp −
r2
1 − α = exp − 2
2σr
log(1 − α) = −

r2
2σr2

r2
2σr2

(A.8)
(A.9)
(A.10)

−2σr2 log(1 − α) = r 2

(A.11)

−2σr2 log(1 − α) = r .

(A.12)

Since α ∈ (0, 1), 1 − α ∈ (0, 1) as well and we can substitute 1 − α with α, giving us

r =

−2σr2 log(α).

(A.13)

From this we can go back to our initial x,y coordinates using x = r cos θ, y = r sin θ and

124

the final result is:

x=

−2σr2 ln α cos(2πφ)

(A.14)

y=

−2σr2 ln α sin(2πφ),

(A.15)

where α and φ are random numbers generated uniformly in the interval (0, 1). This is the
well known Box-Muller transform.

A.3

Uniform laser profile

Now we use the same procedure for the distribution consisting of a uniform circle in the xy
plane, written in radial coordinates as

f (r) =





C =

1
πR02




0

r ≤ R0
.

(A.16)

r > R0

We note that similarly to the Gaussian case, the probability distribution does not depend
on the angle θ, which can then simply be taken a from uniform distribution (0, 2π).
As before, we calculate the cumulative radial probability distribution
2π

g(r ) =

r

dθ
0

r

drrf (r) =
0

dr
0

2r
r2
=
,
R2
R2

(A.17)

and then invert this to get the function that maps the uniform distribution into our desired

125

one
r2
R2

α = g(r ) =
√
r = R α.

In the end the initial x,y coordinates are given by
√
x = R α cos(2πφ)

(A.18)

√
y = R α sin(2πφ),

(A.19)

where α and φ are random numbers generated uniformly in the interval (0, 1).

A.4

Elliptical laser profile

The final case we wish to consider is that of a laser with an elliptical profile, given by the
distribution
f (r) =

3
2πR2

1−

r
R

2

0 ≤ r ≤ R.

(A.20)

Since once again this does not depend on θ, the angular part is easily treated with a uniform
distribution.

126

For the radial part, we integrate f (r) to obtain the cumulative distribution g(r )
2π

g(r ) =

r

drrf (r)

dθ
0

0
r

3
= 2
R
=−

drr

1−

0
2

(R − r2 )3/2
R3
2

=1−

r
R

2

r

0
2 3/2

(R − r )
R3

.

We now follow the usual procedure to invert g(r )

α = g(r ) = 1 −

(R2 − r 2 )3/2
R3

(1 − α)R3 = (R2 − r 2 )3/2
(1 − α)2/3 R2 = R2 − r 2
r 2 = R2 1 − (1 − α)2/3
r =R

1 − (1 − α)2/3 ,

and by using the same substitution used in the Gaussian case, (1 − α) → α, the equations
for the initial coordinates become

x=R

1 − α2/3 cos(2πφ)

(A.21)

y=R

1 − α2/3 sin(2πφ),

(A.22)

where α and φ are random numbers generated uniformly in the interval (0, 1).

127

Appendix B:
Overview of Density Functional
Theory
In this appendix we will give a brief overview of the theory and approximations that govern
Density Functional Theory. For an in depth discussion of both the fundamentals and the
applications of DFT to a variety of systems, the reviews by Gillan [111] and Hafner [112]
offer an solid starting point. Also the Nobel prize talks by Kohn [113] and Pople [114] give
an excellent discussion of the development of methods to simulate the electronic structure
of materials and of the derivation of the DFT equations. More information on the technical
aspects of the computation and its implementation in computational codes can be found in
the review by Cramer and Truhlar [115] and in the papers by Kresse et al. [98–101, 116],
Bloechl [103] and Perdew et al. [102].

B.1

From wavefunctions to a density functional

In principle, if one wishes to study materials using an ab-initio method, all that is required
is to write out the Schroedinger equation for the system and solve it. The resulting ground

128

state energy and wavefunctions uniquely determine all the properties of the system under
consideration. The term ab-initio here is taken to mean methods in which, in addition to the
nuclear mass and coordinates, only the fundamental physical constants enter the calculation
and no empirical parameters are included in the solution.
For a system of Nel electrons and Nat nuclei, the non-relativistic and time independent
Hamiltonian and Schroedinger equations can be written as

ˆ =− 1
H
2mel

Nel

∇2i −
i=1

1
2

Nat

i=1

∇2i
+
Mi

Nel Nel

i=1 i=j

1
+
|ri − rj |

Nat Nat

i=1 j=i

Zi Zj
−
|Ri − Rj |

Nat Nel

i=1 j=1

Zi
|Ri − rj |
(B.23)

ˆ α = Eα Ψα ,
HΨ

(B.24)

where the first two terms are respectively the kinetic energy of the electrons of mass mel
and that of the nuclei of mass Mi . The following terms include electrostatic interactions
between the electrons and nuclei with spatial coordinates ri and Ri . Ψα is the wavefunction
corresponding to the alphath state of the system.
Since we are mostly interested in properties arising from the electronic distribution and
nuclei are approximately 104 times heavier than electrons, with a size that is approximately
10−5 times the size of the electronic orbit, we integrate out the nuclear motion and treat
them as stationary point particles. This is the so-called Born-Oppenheimer approximation,
using which we rewrite the Hamiltonian as

ˆ el = − 1
H
2mel
1
=−
2mel

Nel

Nel Nel

∇2i +
i=1

i=1 i=j

Nel

Nel Nel

∇2i
i=1

+
i=1 i=j

1
−
|ri − rj |
1
−
|ri − rj |

129

Nat Nel

i=1 j=1

Zi
|Ri − rj |

(B.25)

Nel

Vext (ri ),
i=1

(B.26)

and
ˆ el ψel = Eel ψel .
H

(B.27)

The problem we face now is solving this Schroedinger equation to find the antisymmetric
electronic wavefunction ψel , which depends on the positions and spins of the Nel electrons
for fixed nuclear positions Ri . This is a prohibitive task for any system consisting of more
than O(10) electrons [113], leading to the idea of describing the system using the density as
it only depends on 3 coordinates regardless of the size of the system:

n(r) = Nel

···

∗
ψel
(r, r2 , · · · rNel ; σ) ψel (r, r2 , · · · rNel ; σ) dr2 · · · drNel dσ,

(B.28)

where the integrals are over the Nel − 1 spatial degrees of freedom and the spins (indicated
collectively by σ).
Hohenberg and Kohn [117] showed that the ground state density uniquely determines the
Hamiltonian of the system and its ground state energy E0 and that therefore we can rewrite
the problem posed by equations B.26 and B.27 as:

E0 = min E[n(r)]

(B.29)

n(r)

E[n(r)] = Ekin [n] + Eee [n] +

Vext (r)n(r)dr,

(B.30)

where the unknown is now represented by the density n(r).

B.2

The Kohn-Sham equations

Up to this point, the equations written are exact and no approximations have been made.
They are also impossible to solve, as the kinetic and electron-electron contributions to the
130

total energy in Eq. B.30 can’t be written in an exact form as functionals of the density.
To get around this it is necessary to find approximate expressions for Ekin [n] and Eee [n].
Following the work done by Kohn and Sham [118], we start by rewriting the electron-electron
interaction as a sum of the Coulomb interaction term and two non-classical terms that include
the electron exchange and correlation effects:

Eee [n] = ECoul [n] + EX [n] + EC [n]
=

(B.31)

n(r)n(r )
drdr + EX [n] + EC [n].
|r − r |

1
2

(B.32)

We can now rewrite the kinetic energy term as the kinetic energy of independent particles,
Ekin,KS [n], plus the difference of this from the original kinetic energy:

Ekin [n] = Ekin,KS [n] + (Ekin [n] − Ekin,KS [n]).

Putting all of this together gives:

E[n] = Ekin [n] + Eee [n] +
= Ekin,KS [n] +

1
2

Vext (r)n(r)dr
n(r)n(r )
drdr +
|r − r |

Vext (r)n(r)dr

+ EX [n] + EC [n] + Ekin [n] − Ekin,KS [n]
= Ekin,KS [n] +
= Ekin,KS [n] +

1
2

n(r)n(r )
drdr +
|r − r |
1
2

Vext (r)n(r)dr + EXC [n]

n(r )
dr + Vext (r) + VXC (r) n(r)dr,
|r − r |

131

(B.33)

where the last three terms in the second row have been grouped into the unknown functional

EXC = EX [n] + EC [n] + Ekin [n] − Ekin,KS [n] =

VXC (r)n(r)dr.

(B.34)

Using the fact that Ekin,KS [n] corresponds to the kinetic energy of independent particles,
we can now treat the system as being composed of Nel independent particles such that:
Nel

ˆ KS (ri )
h

ˆ ef f =
H

(B.35)

i=1

ˆ KS ψj (r) =
h
VKS =

1
2

1
− ∇2 + VKS (r) ψj (r) =
2

j ψj (r)

n(r )
dr + Vext (r) + VXC (r).
|r − r |

(B.36)
(B.37)

This form of VKS is chosen so that the density obtained from this system of non-interacting
fictitious particles is the same as the ground state density in the true system of interacting
electrons. Approximations now enter our calculations in the choice of the expression for
VXC (r).
A couple of notes to conclude the discussion:
• The Kohn-Sham single particle orbitals, ψj , should not be considered representative
of real electron orbitals in the system. Rather they are a mathematical construction
used to give the correct ground state density. Also the independent particles used in
this description are Kohn-Sham fictitious particles, not real electrons.
• Since VKS depends on the density, it is necessary to do a self-consistent calculation
in which one starts from a guess for the orbitals ψj , then solves the Schroedinger
equation (Eq. B.36) to calculate the new orbitals ψjnew and iterates this procedure

132

until convergence is reached.
• In the work presented in chapter 4, the exchange-correlation functional, EXC , was
described using the generalized gradient approximation (GGA) as implemented by
Perdew, Burke and Ernzerhof (PBE) [102]. A discussion of the various choices that
can be made in terms of exchange-correlation functionals is beyond the scope of this
overview, but can be found in numerous reviews (see for example [115]).
• In addition to the approximation resulting from the choice of EXC or, equivalently,
VXC , the single particle wavefunction is also expanded on a set of basis functions that
must be chosen appropriately. In the work presented here, we used plane waves with
the projector augmented wave method [103] as the basis set.

133

BIBLIOGRAPHY

134

BIBLIOGRAPHY
[1] M. Bovet, S. van Smaalen, H. Berger, R. Gaal, L. Forr´o, L. Schlapbach, and P. Aebi,
“Interplane coupling in the quasi-two-dimensional 1T-TaS2 ,” Phys. Rev. B, vol. 67,
p. 125105, 2003.
[2] C. M. Dobson, “Protein folding and misfolding,” Nature, vol. 426, no. 6968, pp. 884–
890, 2003.
[3] H. Marciniak, M. Fiebig, M. Huth, S. Schiefer, B. Nickel, F. Selmaier, and S. Lochbrunner, “Ultrafast exciton relaxation in microcrystalline pentacene films,” Phys. Rev. Lett.,
vol. 99, p. 176402, Oct 2007.
[4] H. N. Chapman, P. Fromme, A. Barty, T. A. White, R. A. Kirian, A. Aquila, M. S.
Hunter, J. Schulz, D. P. DePonte, U. Weierstall, R. B. Doak, F. R. N. C. Maia, A. V.
Martin, I. Schlichting, L. Lomb, N. Coppola, R. L. Shoeman, S. W. Epp, R. Hartmann, D. Rolles, A. Rudenko, L. Foucar, N. Kimmel, G. Weidenspointner, P. Holl,
M. Liang, M. Barthelmess, C. Caleman, S. Boutet, M. J. Bogan, J. Krzywinski,
C. Bostedt, S. Bajt, L. Gumprecht, B. Rudek, B. Erk, C. Schmidt, A. Homke, C. Reich,
D. Pietschner, L. Struder, G. Hauser, H. Gorke, J. Ullrich, S. Herrmann, G. Schaller,
F. Schopper, H. Soltau, K.-U. Kuhnel, M. Messerschmidt, J. D. Bozek, S. P. HauRiege, M. Frank, C. Y. Hampton, R. G. Sierra, D. Starodub, G. J. Williams, J. Hajdu, N. Timneanu, M. M. Seibert, J. Andreasson, A. Rocker, O. Jonsson, M. Svenda,
S. Stern, K. Nass, R. Andritschke, C.-D. Schroter, F. Krasniqi, M. Bott, K. E. Schmidt,
X. Wang, I. Grotjohann, J. M. Holton, T. R. M. Barends, R. Neutze, S. Marchesini,
R. Fromme, S. Schorb, D. Rupp, M. Adolph, T. Gorkhover, I. Andersson, H. Hirsemann, G. Potdevin, H. Graafsma, B. Nilsson, and J. C. H. Spence, “Femtosecond x-ray
protein nanocrystallography,” Nature, vol. 470, pp. 73–77, Feb. 2011.
[5] A. H. Zewail, “Four-dimensional electron microscopy,” Science, vol. 328, no. 5975,
pp. 187–193, 2010.
[6] C.-Y. Ruan, F. Vigliotti, V. A. Lobastov, S. Chen, and A. H. Zewail, “Ultrafast electron
crystallography: Transient structures of molecules, surfaces, and phase transitions,”
Proc. Natl. Acad. Sci. U.S.A., vol. 101, no. 5, pp. 1123–1128, 2004.
[7] A. Lindenberg, J. Larsson, K. Sokolowski-Tinten, K. Gaffney, C. Blome, O. Synnergren, J. Sheppard, C. Caleman, A. MacPhee, D. Weinstein, D. Lowney, T. Allison,
T. Matthews, R. Falcone, A. Cavalieri, D. Fritz, S. Lee, P. Bucksbaum, D. Reis,
135

J. Rudati, P. Fuoss, C. Kao, D. Siddons, R. Pahl, J. Als-Nielsen, S. Duesterer, R. Ischebeck, H. Schlarb, H. Schulte-Schrepping, T. Tschentscher, J. Schneider, D. von der
Linde, O. Hignette, F. Sette, H. Chapman, R. Lee, T. Hansen, S. Techert, J. Wark,
M. Bergh, G. Huldt, D. van der Spoel, N. Timneanu, J. Hajdu, R. Akre, E. Bong,
P. Krejcik, J. Arthur, S. Brennan, K. Luening, and J. Hastings, “Atomic-scale visualization of inertial dynamics,” Science, vol. 308, no. 5720, pp. 392–395, 2005.
[8] P. Baum, D. S. Yang, and A. H. Zewail, “4d visualization of transitional structures in
phase transformations by electron diffraction,” Science, vol. 318, no. 5851, pp. 788–792,
2007.
[9] M. Harb, R. Ernstorfer, C. T. Hebeisen, G. Sciaini, W. Peng, T. Dartigalongue, M. A.
Eriksson, M. G. Lagally, S. G. Kruglik, and R. J. D. Miller, “Electronically driven
structure changes of si captured by femtosecond electron diffraction,” Phys. Rev. Lett.,
vol. 100, no. 15, 2008.
[10] R. K. Raman, R. A. Murdick, R. J. Worhatch, Y. Murooka, S. D. Mahanti, T. R. T.
Han, and C. Y. Ruan, “Electronically driven fragmentation of silver nanocrystals
revealed by ultrafast electron crystallography,” Phys. Rev. Lett., vol. 104, no. 12,
p. 123401, 2010.
[11] P. Musumeci, J. T. Moody, C. M. Scoby, M. S. Gutierrez, and M. Westfall, “Laserinduced melting of a single crystal gold sample by time-resolved ultrafast relativistic
electron diffraction,” Appl. Phys. Lett., vol. 97, no. 6, 2010.
[12] R. A. Murdick, R. K. Raman, Y. Murooka, and C. Y. Ruan, “Photovoltage dynamics
of the hydroxylated si(111) surface investigated by ultrafast electron diffraction,” Phys.
Rev. B, vol. 77, no. 24, p. 245329, 2008.
[13] C. Y. Ruan, Y. Murooka, R. K. Raman, R. A. Murdick, R. J. Worhatch, and A. Pell,
“The development and applications of ultrafast electron nanocrystallography,” Microsc. Microanal., vol. 15, no. 4, pp. 323–337, 2009.
[14] R. D. Averitt and A. Taylor, “Ultrafast optical and far-infrared quasiparticle dynamics in correlated electron materials,” J. Phys.: Condens. Matter, vol. 14, no. 50,
pp. R1357–R1390, 2002.
[15] Z. Yang, C. Ko, and S. Ramanathan, “Oxide electronics utilizing ultrafast metalinsulator transitions,” Annu. Rev. Mater. Res., vol. 41, no. 1, pp. 337–367, 2011.

136

[16] Z. Tao, T.-R. T. Han, S. D. Mahanti, P. M. Duxbury, F. Yuan, C.-Y. Ruan, K. Wang,
and J. Wu, “Decoupling of structural and electronic phase transitions in VO2 ,” Phys.
Rev. Lett., vol. 109, p. 166406, 2012.
[17] A. H. Zewail, “Femtochemistry: Atomic-scale dynamics of the chemical bond,” J.
Phys. Chem. A, vol. 104, no. 24, pp. 5660–5694, 2000.
[18] S. Bratos, J. C. Leicknam, S. Pommeret, and G. Gallot, “Laser spectroscopic visualization of hydrogen bond motions in liquid water,” J. Mol. Struct., vol. 708, no. 1-3,
pp. 197–203, 2004.
[19] C. Y. Ruan, V. A. Lobastov, F. Vigliotti, S. Y. Chen, and A. H. Zewail, “Ultrafast
electron crystallography of interfacial water,” Science, vol. 304, no. 5667, pp. 80–84,
2004.
[20] R. A. Nicodemus, S. A. Corcelli, J. L. Skinner, and A. Tokmakoff, “Collective hydrogen
bond reorganization in water studied with temperature-dependent ultrafast infrared
spectroscopy,” J. Phys. Chem. B, vol. 115, no. 18, pp. 5604–5616, 2011.
[21] M. D. Feit, A. M. Komashko, and A. M. Rubenchik, “Ultra-short pulse laser interaction
with transparent dielectrics,” Appl Phys A-mater, vol. 79, no. 7, pp. 1657–1661, 2004.
[22] M. Murillo, “Ultrafast dynamics of strongly coupled plasmas,” Phys. Rev. Lett., vol. 96,
no. 16, 2006.
[23] B. Dromey, S. Rykovanov, M. Yeung, R. Hoerlein, D. Jung, D. C. Gautier, T. Dzelzainis, D. Kiefer, S. Palaniyppan, R. Shah, J. Schreiber, H. Ruhl, J. C. Fernandez, C. L. S.
Lewis, M. Zepf, and B. M. Hegelich, “Coherent synchrotron emission from electron
nanobunches formed in relativistic laser-plasma interactions,” Nat. Phys., vol. 8, no. 11,
pp. 804–808, 2012.
[24] S. P. Hau-Riege, A. Graf, T. Doeppner, R. A. London, J. Krzywinski, C. Fortmann, S. H. Glenzer, M. Frank, K. Sokolowski-Tinten, M. Messerschmidt, C. Bostedt,
S. Schorb, J. A. Bradley, A. Lutman, D. Rolles, A. Rudenko, and B. Rudek, “Ultrafast transitions from solid to liquid and plasma states of graphite induced by x-ray
free-electron laser pulses,” Phys. Rev. Lett., vol. 108, no. 21, 2012.
[25] W. E. King, G. H. Campbell, A. Frank, B. Reed, J. F. Schmerge, B. J. Siwick, B. C.
Stuart, and P. M. Weber, “Ultrafast electron microscopy in materials science, biology,
and chemistry,” J. Appl. Phys., vol. 97, no. 11, p. 111101, 2005.

137

[26] J. Portman, H. Zhang, Z. Tao, K. Makino, M. Berz, P. M. Duxbury, and C.-Y.
Ruan, “Computational and experimental characterization of high-brightness beams for
femtosecond electron imaging and spectroscopy,” Appl. Phys. Lett., vol. 103, no. 25,
p. 253115, 2013.
[27] S. Hellmann, T. Rohwer, M. Kallne, K. Hanff, C. Sohrt, A. Stange, A. Carr, M. Murnane, H. Kapteyn, and L. Kipp, “Time-domain classification of charge-density-wave
insulators,” Nat. Commun., p. 1069, 2012.
[28] B. Sipos, A. F. Kusmartseva, A. Akrap, H. Berger, L. Forro, and E. Tutis, “From mott
state to superconductivity in 1T-TaS2 ,” Nat. Mater., vol. 7, no. 12, pp. 960–965, 2008.
[29] W. Setyawan and S. Curtarolo, “High-throughput electronic band structure calculations: Challenges and tools,” Comput. Mater. Sci., vol. 49, no. 2, pp. 299 – 312, 2010.
[30] D. J. Flannigan and A. H. Zewail, “4d electron microscopy: Principles and applications,” Acc. Chem. Res., vol. 45, no. 10, pp. 1828–1839, 2012.
[31] T. Brabec and F. Krausz, “Intense few-cycle laser fields: Frontiers of nonlinear optics,”
Rev. Mod. Phys., vol. 72, pp. 545–591, Apr 2000.
[32] G. A. Mourou, T. Tajima, and S. V. Bulanov, “Optics in the relativistic regime,” Rev.
Mod. Phys., vol. 78, pp. 309–371, Apr 2006.
[33] P. B. Corkum and F. Krausz, “Attosecond science,” Nat. Phys., vol. 3, pp. 381–387,
June 2007.
[34] M. Chergui and A. H. Zewail, “Electron and x-ray methods of ultrafast structural
dynamics: Advances and applications,” Chem. Phys. Chem., vol. 10, no. 1, pp. 28–43,
2009.
[35] F. Carbone, P. Musumeci, O. Luiten, and C. Hebert, “A perspective on novel sources
of ultrashort electron and x-ray pulses,” Chemical Physics, vol. 392, no. 1, pp. 1–9,
2012.
[36] Zhang, The Fast Multipole Method In The Differential Algebra Framework For The
Calculation Of 3d Space Charge Fields. PhD thesis, Michigan State University, 2012.
[37] A. M. Michalik and J. E. Sipe, “Analytic model of electron pulse propagation in ultrafast electron diffraction experiments,” J. Appl. Phys., vol. 99, no. 5, p. 054908,
2006.
138

[38] A. M. Michalik and J. E. Sipe, “Evolution of non-gaussian electron bunches in ultrafast electron diffraction experiments: Comparison to analytic model,” J. Appl. Phys.,
vol. 105, no. 8, p. 084913, 2009.
[39] J. A. Berger and W. A. Schroeder, “Semianalytic model of electron pulse propagation:
Magnetic lenses and rf pulse compression cavities,” J. Appl. Phys., vol. 108, no. 12,
p. 124905, 2010.
[40] P. Hawkes and J. Spence, Science of Microscopy. Springer, 2008.
[41] M. Aidelsburger, F. O. Kirchner, F. Krausz, and P. Baum, “Single-electron pulses for
ultrafast diffraction,” Proc. Natl. Acad. Sci. U.S.A., vol. 107, no. 46, pp. 19714–19719,
2010.
[42] D. S. Yang, O. F. Mohammed, and A. H. Zewail, “Scanning ultrafast electron microscopy,” Proc. Natl. Acad. Sci. U.S.A., vol. 107, no. 34, pp. 14993–14998, 2010.
[43] A. Zewail and J. Thomas, 4D Electron Microscopy: Imaging in Space and Time. Imperial College Press, 2010.
[44] G. H. Kassier, N. Erasmus, K. Haupt, I. Boshoff, R. Siegmund, S. M. M. Coelho, and
H. Schwoerer, “Photo-triggered pulsed cavity compressor for bright electron bunches
in ultrafast electron diffraction,” Appl. Phys. B, vol. 109, no. 2, pp. 249–257, 2012.
[45] G. F. Mancini, B. Mansart, S. Pagano, B. van der Geer, M. de Loos, and F. Carbone,
“Design and implementation of a flexible beamline for fs electron diffraction experiments,” Nucl. Instrum. Methods Phys. Res., Sect. A, vol. 691, pp. 113–122, 2012.
[46] T. van Oudheusden, E. F. de Jong, S. B. van der Geer, W. P. E. M. O. Root, O. J.
Luiten, and B. J. Siwick, “Electron source concept for single-shot sub-100 fs electron
diffraction in the 100 kev range,” J. Appl. Phys., vol. 102, no. 9, p. 093501, 2007.
[47] I. V. Bazarov, B. M. Dunham, and C. K. Sinclair, “Maximum achievable beam brightness from photoinjectors,” Phys. Rev. Lett., vol. 102, p. 104801, 2009.
[48] D. H. Dowell and J. F. Schmerge, “Quantum efficiency and thermal emittance of metal
photocathodes,” Phys. Rev. Spec. Top. Accel Beams, vol. 12, no. 7, p. 074201, 2009.
[49] D. H. Dowell, I. Bazarov, B. Dunham, K. Harkay, C. Hernandez-Garcia, R. Legg,
H. Padmore, T. Rao, J. Smedley, and W. Wan, “Cathode r&d for future light sources,”
Nucl. Instrum. Methods Phys. Res., Sect. A, vol. 622, no. 3, pp. 685–697, 2010.
139

[50] K. Floettmann, “Some basic features of the beam emittance,” Phys. Rev. Spec. Top.
Accel Beams, vol. 6, no. 3, p. 034202, 2003.
[51] K. L. Jensen, N. A. Moody, D. W. Feldman, E. J. Montgomery, and P. G. O’Shea,
“Photoemission from metals and cesiated surfaces,” J. Appl. Phys., vol. 102, no. 7,
p. 074902, 2007.
[52] K. L. Jensen, P. G. O’Shea, D. W. Feldman, and J. L. Shaw, “Emittance of a field
emission electron source,” J. Appl. Phys., vol. 107, no. 1, p. 014903, 2010.
[53] O. Luiten, S. van der Geer, M. de Loos, F. Kiewiet, and M. van der Wiel, “How
to realize uniform three-dimensional ellipsoidal electron bunches,” Phys. Rev. Lett.,
vol. 93, no. 9, 2004.
[54] K. Makino and M. Berz, “Cosy infinity version 9,” Nucl. Instrum. Methods Phys. Res.,
Sect. A, vol. 558, pp. 346–350, 2006.
[55] C. N. Berglund and W. E. Spicer, “Photoemission studies of copper and silver: Theory,” Phys. Rev., vol. 136, pp. A1030–A1044, 1964.
[56] T. Srinivasan-Rao, J. Fischer, and T. Tsang, “Photoemission studies on metals using
picosecond ultraviolet laser pulses,” J. Appl. Phys., vol. 69, no. 5, pp. 3291–3296, 1991.
[57] X. Jiang, C. N. Berglund, A. E. Bell, and W. A. Mackie, “Os3: Photoemission from
gold thin films for application in multiphotocathode arrays for electron beam lithography,” Papers from the 42nd international conference on electron, ion, and photon
beam technology and nanofabrication, vol. 16, no. 6, pp. 3374–3379, 1998.
[58] L. Rangarajan and G. Bhide, “Photoemission energy distribution studies of gold thin
films under uv excitation by a photoelectron spectroscopic method,” Vacuum, vol. 30,
pp. 515 – 522, 1980.
[59] R. Srinivasan, V. Lobastov, C.-Y. Ruan, and A. Zewail, “Ultrafast electron diffraction
(ued),” Helv. Chim. Acta, vol. 86, no. 6, pp. 1761–1799, 2003.
[60] R. Akre, D. Dowell, P. Emma, J. Frisch, S. Gilevich, G. Hays, P. Hering, R. Iverson, C. Limborg-Deprey, H. Loos, A. Miahnahri, J. Schmerge, J. Turner, J. Welch,
W. White, and J. Wu, “Commissioning the linac coherent light source injector,” Phys.
Rev. ST Accel. Beams, vol. 11, p. 030703, Mar 2008.

140

[61] H. Zhang and M. Berz, “The fast multipole method in the differential algebra framework,” Nucl. Instrum. Methods Phys. Res., Sect. A, vol. 645, no. 1, pp. 338 – 344,
2011.
[62] J. Jackson, Classical electrodynamics. Wiley, 1975.
[63] W. Schroeder, J. Berger, A. Nicholls, N. Zaluzec, J. Hiller, and D. Miller, “High brightness nano-patterned photocathode electron sources for uem,” Microsc. Microanal.,
vol. 16, pp. 492–493, 7 2010.
[64] P. Musumeci, J. T. Moody, R. J. England, J. B. Rosenzweig, and T. Tran, “Experimental generation and characterization of uniformly filled ellipsoidal electron-beam
distributions,” Phys. Rev. Lett., vol. 100, no. 24, 2008.
[65] T. van Oudheusden, P. L. E. M. Pasmans, S. B. van der Geer, M. J. de Loos,
M. J. van der Wiel, and O. J. Luiten, “Compression of subrelativistic space-chargedominated electron bunches for single-shot femtosecond electron diffraction,” Phys.
Rev. Lett., vol. 105, no. 26, 2010.
[66] Z. Tao, H. Zhang, P. M. Duxbury, M. Berz, and C.-Y. Ruan, “Space charge effects in
ultrafast electron diffraction and imaging,” J. Appl. Phys., vol. 111, no. 4, p. 044316,
2012.
[67] A. Valfells, D. W. Feldman, M. Virgo, P. G. O’Shea, and Y. Y. Lau, “Effects of pulselength and emitter area on virtual cathode formation in electron guns,” Phys. Plasmas,
vol. 9, no. 5, pp. 2377–2382, 2002.
[68] N. Ravishanker and D. Dey, A First Course in Linear Model Theory. Chapman &
Hall/CRC Texts in Statistical Science, Taylor & Francis, 2001.
[69] P. Pasmans, G. van den Ham, S. D. Conte, S. van der Geer, and O. Luiten, “Microwave
tm010 cavities as versatile 4d electron optical elements,” Ultramicroscopy, vol. 127,
no. 0, pp. 19 – 24, 2013. Frontiers of Electron Microscopy in Materials Science.
[70] L. Landau, The Classical Theory of Fields. Course Of Theoretical Physics, Elsevier
Science, 1980.
[71] T. Dray, The Geometry of Special Relativity. Taylor & Francis, 2012.
[72] T. S. V. W. F. B. Press, W.H., Numerical Recipes 3rd Edition: The Art of Scientific
Computing. Cambridge University Press, 2007.
141

[73] B. K. Vollhardt Dieter and K. Marcus, Strongly Correlated Systems: Theoretical Methods, vol. 171 of Springer Series in Solid-State Sciences. Springer Berlin Heidelberg,
2012. ISBN 978-3-642-21830-9.
[74] R. Peierls, Quantum Theory of Solids. International Series of Monographs on Physics,
Clarendon Press, 1955.
[75] G. Gruner, Density Waves In Solids. Frontiers in Physics Series, Addison-Wesley,
1994.
[76] A. R. Beal, H. P. Hughes, and W. Y. Liang, “The reflectivity spectra of some group va
transition metal dichalcogenides,” J. Phys. C: Solid State Phys., vol. 8, no. 24, p. 4236,
1975.
[77] R. E. Thomson, B. Burk, A. Zettl, and J. Clarke, “Scanning tunneling microscopy of
the charge-density-wave structure in 1T-TaS2 ,” Phys. Rev. B, vol. 49, pp. 16899–16916,
1994.
[78] K. Fung, J. Steeds, and J. Eades, “Application of convergent beam electron diffraction
to study the stacking of layers in transition-metal dichalcogenides,” Physica B+C,
vol. 99, pp. 47 – 50, 1980.
[79] S. Tanda, T. Sambongi, T. Tani, and S. Tanaka, “X-ray study of charge density wave
structure in 1T-TaS2 ,” J. Phys. Soc. Jpn., vol. 53, no. 2, pp. 476–479, 1984.
[80] K. Nakanishi and H. Shiba, “Theory of three-dimensional orderings of charge-density
waves in 1T-TaX2 (X: S, Se),” J. Phys. Soc. Jpn., vol. 53, no. 3, pp. 1103–1113, 1984.
[81] N. V. Smith, S. D. Kevan, and F. J. DiSalvo, “Band structures of the layer compounds
1T-TaS2 and 2H-TaSe2 in the presence of commensurate charge-density waves,” J.
Phys. C: Solid State Phys., vol. 18, no. 16, p. 3175, 1985.
[82] T. Pillo, J. Hayoz, H. Berger, M. Grioni, L. Schlapbach, and P. Aebi, “Remnant fermi
surface in the presence of an underlying instability in layered 1T-TaS2 ,” Phys. Rev.
Lett., vol. 83, pp. 3494–3497, Oct 1999.
[83] T. Pillo, J. Hayoz, H. Berger, R. Fasel, L. Schlapbach, and P. Aebi, “Interplay between
electron-electron interaction and electron-phonon coupling near the fermi surface of
1T-TaS2 ,” Phys. Rev. B, vol. 62, pp. 4277–4287, 2000.

142

[84] L. Perfetti, T. A. Gloor, F. Mila, H. Berger, and M. Grioni, “Unexpected periodicity in the quasi-two-dimensional mott insulator 1T-TaS2 revealed by angle-resolved
photoemission,” Phys. Rev. B, vol. 71, p. 153101, Apr 2005.
[85] F. Clerc, C. Battaglia, M. Bovet, L. Despont, C. Monney, H. Cercellier, M. G. Garnier, P. Aebi, H. Berger, and L. Forr´o, “Lattice-distortion-enhanced electron-phonon
coupling and fermi surface nesting in 1T-TaS2 ,” Phys. Rev. B, vol. 74, p. 155114, 2006.
[86] S. Hellmann, M. Beye, C. Sohrt, T. Rohwer, F. Sorgenfrei, H. Redlin, M. Kall¨ane,
M. Marczynski-B¨
uhlow, F. Hennies, M. Bauer, A. F¨ohlisch, L. Kipp, W. Wurth, and
K. Rossnagel, “Ultrafast melting of a charge-density wave in the mott insulator 1TTaS2 ,” Phys. Rev. Lett., vol. 105, p. 187401, Oct 2010.
[87] B. Dardel, M. Grioni, D. Malterre, P. Weibel, Y. Baer, and F. L´evy, “Temperaturedependent pseudogap and electron localization in 1T-TaS2 ,” Phys. Rev. B, vol. 45,
pp. 1462–1465, Jan 1992.
[88] K. Rossnagel and N. V. Smith, “Spin-orbit coupling in the band structure of reconstructed 1T-TaS2 ,” Phys. Rev. B, vol. 73, p. 073106, 2006.
[89] P. Darancet, A. Millis, and C. A. Marianetti, “Three dimensional metallic and two
dimensional insulating behavior in tantalum dichalcogenides,” 2014.
[90] A. Spijkerman, J. L. de Boer, A. Meetsma, G. A. Wiegers, and S. van Smaalen, “X-ray
crystal-structure refinement of the nearly commensurate phase of 1T-TaS2 in (3+2)dimensional superspace,” Phys. Rev. B, vol. 56, pp. 13757–13767, 1997.
[91] F. Zwick, H. Berger, I. Vobornik, G. Margaritondo, L. Forr´o, C. Beeli, M. Onellion,
G. Panaccione, A. Taleb-Ibrahimi, and M. Grioni, “Spectral consequences of broken
phase coherence in 1T-TaS2 ,” Phys. Rev. Lett., vol. 81, pp. 1058–1061, Aug 1998.
[92] A. Y. Liu, “Electron-phonon coupling in compressed 1T-TaS2 : Stability and superconductivity from first principles,” Phys. Rev. B, vol. 79, p. 220515, 2009.
[93] T. Ritschel, J. Trinckauf, G. Garbarino, M. Hanfland, M. v. Zimmermann, H. Berger,
B. B¨
uchner, and J. Geck, “Pressure dependence of the charge density wave in 1T-TaS2
and its relation to superconductivity,” Phys. Rev. B, vol. 87, p. 125135, 2013.
[94] R. Ang, Y. Tanaka, E. Ieki, K. Nakayama, T. Sato, L. J. Li, W. J. Lu, Y. P. Sun, and
T. Takahashi, “Real-space coexistence of the melted mott state and superconductivity
in fe-substituted 1T-TaS2 ,” Phys. Rev. Lett., vol. 109, p. 176403, 2012.
143

[95] L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, H. Berger, S. Biermann,
P. S. Cornaglia, A. Georges, and M. Wolf, “Time evolution of the electronic structure
of 1T-TaS2 through the insulator-metal transition,” Phys. Rev. Lett., vol. 97, p. 067402,
2006.
[96] L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, M. Wolf, H. Berger, S. Biermann, and A. Georges, “Femtosecond dynamics of electronic states in the mott insulator 1T-TaS2 by time resolved photoelectron spectroscopy,” New J. Phys., vol. 10,
no. 5, p. 053019, 2008.
[97] J. K. Freericks, H. R. Krishnamurthy, Y. Ge, A. Y. Liu, and T. Pruschke, “Theoretical description of time-resolved pump/probe photoemission in TaS2 : a single-band
DFT+DMFT(NRG) study within the quasiequilibrium approximation,” Phys. Status
Solidi B, vol. 246, no. 5, pp. 948–954, 2009.
[98] G. Kresse and J. Hafner, “Ab initio molecular dynamics for liquid metals,” Phys. Rev.
B, vol. 47, pp. 558–561, 1993.
[99] G. Kresse and J. Hafner, “Ab initio molecular-dynamics simulation of the liquid-metalamorphous-semiconductor transition in germanium,” Phys. Rev. B, vol. 49, pp. 14251–
14269, 1994.
[100] G. Kresse and J. Furthmueller, “Efficiency of ab-initio total energy calculations for
metals and semiconductors using a plane-wave basis set,” Comput. Mater. Sci., vol. 6,
no. 1, pp. 15 – 50, 1996.
[101] G. Kresse and J. Furthmueller, “Efficient iterative schemes for ab initio total-energy
calculations using a plane-wave basis set,” Phys. Rev. B, vol. 54, pp. 11169–11186,
1996.
[102] J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made
simple,” Phys. Rev. Lett., vol. 77, pp. 3865–3868, 1996.
[103] P. E. Bloechl, “Projector augmented-wave method,” Phys. Rev. B, vol. 50, pp. 17953–
17979, 1994.
[104] P. Pulay, “Convergence acceleration of iterative sequences. the case of scf iteration,”
Chem. Phys. Lett., vol. 73, no. 2, pp. 393 – 398, 1980.
[105] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton,
“Electron-energy-loss spectra and the structural stability of nickel oxide: An lsda+u
study,” Phys. Rev. B, vol. 57, pp. 1505–1509, Jan 1998.
144

[106] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, “Phonons and related crystal properties from density-functional perturbation theory,” Rev. Mod. Phys., vol. 73,
pp. 515–562, Jul 2001.
[107] M. Gajdos, K. Hummer, G. Kresse, J. Furthm¨
uller, and F. Bechstedt, “Linear optical properties in the projector-augmented wave methodology,” Phys. Rev. B, vol. 73,
p. 045112, Jan 2006.
[108] L. Perfetti, A. Georges, S. Florens, S. Biermann, S. Mitrovic, H. Berger, Y. Tomm,
H. H¨ochst, and M. Grioni, “Spectroscopic signatures of a bandwidth-controlled mott
transition at the surface of 1T-TaSe2 ,” Phys. Rev. Lett., vol. 90, p. 166401, Apr 2003.
[109] J. K. Freericks, “Quenching bloch oscillations in a strongly correlated material:
Nonequilibrium dynamical mean-field theory,” Phys. Rev. B, vol. 77, p. 075109, 2008.
[110] M. A. Marques, N. T. Maitra, F. M. Nogueira, E. Gross, and A. Rubio, eds., Fundamentals of Time-Dependent Density Functional Theory, vol. 837 of Lecture Notes in
Physics. Springer Berlin Heidelberg, 2012.
[111] M. J. Gillan, “The virtual matter laboratory,” Contemp. Phys., vol. 38, no. 2, pp. 115–
130, 1997.
[112] J. Hafner, “Atomic-scale computational materials science,” Acta Mater., vol. 48, no. 1,
pp. 71 – 92, 2000.
[113] W. Kohn, “Nobel lecture: Electronic structure of matter - wave functions and density
functionals,” Rev. Mod. Phys., vol. 71, pp. 1253–1266, 1999.
[114] J. A. Pople, “Nobel lecture: Quantum chemical models,” Rev. Mod. Phys., vol. 71,
pp. 1267–1274, 1999.
[115] C. J. Cramer and D. G. Truhlar, “Density functional theory for transition metals and
transition metal chemistry,” Phys. Chem. Chem. Phys., vol. 11, pp. 10757–10816, 2009.
[116] G. Kresse and D. Joubert, “From ultrasoft pseudopotentials to the projector
augmented-wave method,” Phys. Rev. B, vol. 59, pp. 1758–1775, 1999.
[117] P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev., vol. 136,
pp. B864–B871, 1964.

145

[118] W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation
effects,” Phys. Rev., vol. 140, pp. A1133–A1138, 1965.

146