ANALYTICAL AND COMPUTATIONAL STUDY OF ELECTRON

BUNCH DYNAMICS

By

Xukun Xiang

A DISSERTATION

Submitted to

Michigan State University

in partial fulﬁllment of the requirements

for the degree of

Physics – Doctor of Philosophy

2020

ABSTRACT

ANALYTICAL AND COMPUTATIONAL STUDY OF ELECTRON BUNCH

DYNAMICS

By

Xukun Xiang

This dissertation is centered on the analytical and numerical study of ellipsoidal elec-

tron bunch dynamics. We are particularly interested in the focusing process of the probing

electron bunches in the ultrafast electron diﬀraction/microscopy system, so that we can

improve the temporal and spectral resolution of the ultrafast experiments.

More speciﬁcally, to understand the collective space charge eﬀects throughout the

bunch evolution, we employed several analytic models to describe the bunch dynamics.

We start with an extension of the mean-ﬁeld model using ordinary diﬀerential equations.

Analysis of this mean-ﬁeld model leads to the identiﬁcation of a longitudinal critical chirp,

which separates two regimes for particle trajectories for the longitudinal focusing of the

bunches: (1) bounce-back, where the particles reverse their direction at the waist of the

focusing process, and (2) cross-over, where the bunch experiences a singularity where the

bunch width reduces to zero. We show that time can be scaled by the initial plasma

frequency, and that the critical chirp becoming dimensionless and depends solely on the

initial bunch aspect ratio.

In order to study the emittance eﬀect on the bunch dynamics, we introduce the self-

similar analytical (SSA) model, a statistical method describing the second order moments

dynamics to model the evolution of an ellipsoidal electron bunch. We also discussed

its linear chirp assumption, explaining how it is the key assumption that leads to the

emittance conservation according to the SSA model. We discuss the statistical nature of

bunch emittance noting that the space charge eﬀect of the uniform density proﬁle and

of the Gaussian proﬁle are close to each other in the SSA model. The impact from a

changing emittance is captured by an additional term in the modiﬁed SSA model, which

is then equivalent to the Kapchinsky-Vladimirsky (K-V) envelope equation in accelerator

physics. We point out that the application of the statistical methods can extend beyond

the uniform ellipsoidal bunch, while the accuracy of the SSA prediction is mainly related

to the discrepancy between the actual density proﬁle and the uniform density proﬁle.

We present the Molecular Dynamics (MD) simulation results for the longitudinal fo-

cusing process of uniform spheroidal electron bunches. The comparison of the longitudinal

width evolution between the MD simulations and the SSA predictions shows the impact of

a varying emittance on bunch evolution. We propose two competing mechanisms for the

change of emittance throughout the compression process. The disorder-induced heating

(DIH) increases the emittance in both degrees of freedom while the diﬀerence in the SSA

temperature generates emittance transfer between degrees of freedom. In addition, the

non-uniform density proﬁle at the focal point introduces signiﬁcant emittance growth in

both the longitudinal and transverse directions.

Copyright by
XUKUN XIANG
2020

To my dear family

v

ACKNOWLEDGMENTS

Doctoral study is an incredible life experience that is so priceless with the help and

support of many great people. I want to say a massive THANK YOU to all the wonderful

people that I have encountered and learned from during this memorable journey.

First, I would like to express my greatest gratitude to my advisor, Dr. Phil Duxbury,

for all his guidance, support, and tolerance. He taught me how to scientiﬁcally solve

a problem as a properly trained scientist and researcher. I would like to thank all my

committee members for the great advice and guidance for improving my research. I would

like to thank Dr. Steven Lund and Dr. Chong-Yu Ruan. I really appreciate the experience

of our collaboration and the opportunity to learn from you. I also thank Dr. Reginald

Ronningen, Dr. David Tomanek, and Dr. Carl Schmidt for taking the time to understand

my research. My Ph.D. work cannot be done without collaboration with my colleague,

Dr. Brandon Zerbe, for his generous help and inspiring discussions.

I would like to extend my sincere gratitude to my friends along this journey. I am so

lucky to share wonderful memories and delicious food with Faran Zhou, Zhensheng Tao,

Mengze Zhu, Jie Guan, Xu Lu, Didi Luo, Xueying Huyan, Nan Du, Zhiyi Su and many

more. I would also like to send special thanks to Dr. Ken Frank, Dr. Andy Anderson,

and Zen Zhong for all the inspiring discussions that expand my understanding of life and

myself.

I would like to thank my parents for their unconditional love and support.

I am

extremely grateful that you stopped asking me how much longer it will take me to ﬁnish

after my ﬁfth year. I love you!

vi

Finally, I would like to give my most special thanks to my wife, Dr. Qinyun Lin, who

ﬁnished her Ph.D. a lot faster than me. Thank you for always being there with me. I am

the luckiest person in the world to share my life with you!

vii

TABLE OF CONTENTS

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

x

KEY TO ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii

Chapter 1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 The Ultrafast Electron Microscopy . . . . . . . . . . . . . . . . . . . . . .
1.2
Importance of our analytical modeling . . . . . . . . . . . . . . . . . . . .
1.3 The Existing analytical models
. . . . . . . . . . . . . . . . . . . . . . . .
1.4 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
4
5
7

2.1 Mean-ﬁeld Theory for a uniform spheroidal bunch . . . . . . . . . . . . . .

Chapter 2 Mean-ﬁeld theory for dynamics of ellipsoidal electron bunches 9
9
2.1.1 Mean-ﬁeld theory for a uniform ellipsoidal bunch . . . . . . . . . . 10
2.1.2
Self-Similar Evolution and Equations of Motion . . . . . . . . . . . 16
Initial condition and Coulomb Explosion . . . . . . . . . . . . . . . 18
2.1.3
2.1.4 Critical Chirp ω∗
c , Crossover and Bounce-back . . . . . . . . . . . . 19
2.2 Coulomb Potential inside a Gaussian Bunch . . . . . . . . . . . . . . . . . 25
2.2.1 Potential Contribution from one single shell
. . . . . . . . . . . . . 26
2.2.2 Potential at an interior point of a Gaussian bunch . . . . . . . . . . 29
2.3 Discussion about the linearity assumption of MFT . . . . . . . . . . . . . . 32

Chapter 3

Statistical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.1 The Self-similar Analytical Model
. . . . . . . . . . . . . . . . . . . . . . . 34
3.2 Bunch evolution with conserved emittance . . . . . . . . . . . . . . . . . . 40
3.3 Linear chirp assumption conserves emittance . . . . . . . . . . . . . . . . . 43
3.3.1 The Linear chirp assumption in detail . . . . . . . . . . . . . . . . . 44
3.3.2 A Comparison Between the SSA model and the Analytic Gaussian

model

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.3.3 Discussion about AG . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
linear chirp and self-similar evolution assumptions . . . . . . . . . . 59
3.4.1
3.4.2 Kinetic energy transfer between degrees of freedom . . . . . . . . . 60
Statistical discussion and the K-V envelope equations . . . . . . . . 61
3.4.3
3.4.4
SSA with emittance evolution . . . . . . . . . . . . . . . . . . . . . 65

Chapter 4 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . 67
4.1 Simulation method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.2 Warm initial condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.3 Width evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.4 Emittance evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

viii

4.5.1 Density proﬁle prediction from SSA and its linear chirp assumption

4.4.1 Disorder-induced heating (DIH) . . . . . . . . . . . . . . . . . . . . 77
4.4.2 Emittance transfer between degrees of freedom . . . . . . . . . . . . 78
4.4.3 Phenomenological description of the emittance evolution . . . . . . 82
4.5 Density proﬁle at tc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
83
4.6 Bunch phase space display for cylindrical bunch . . . . . . . . . . . . . . . 87
4.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.7.1 More on the emittance evolution . . . . . . . . . . . . . . . . . . . 89
4.7.2
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.7.3 Comparing with the experiment values . . . . . . . . . . . . . . . . 91

Sampling error

Chapter 5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

ix

LIST OF FIGURES

Figure 1.1: Ultrafast electron beam column with optics for controlling electron
bunch phase space evolution. (a) Schematic drawing of a prototype
UEM beamline. (b)–(d) The phase space evolution in the injector por-
tion of the beamline (enclosed in dashed line). (e) The conceptual out-
line of the atomic grating approach to characterize the energy spread of
the electron bunches. Adapted from “Active control of bright electron
beams with RF optics for femtosecond microscopy” by Williams, J. et.
al., 2017, Structural Dynamics, 4(4), 044035. . . . . . . . . . . . . . . .

3

Figure 2.1: Longitudinal width evolution Z(τ ) of prolate spheroids with (α0 =
10/75) driven by diﬀerent initial chirps in numeric solutions of MFT
c ) to well above (−2.0ω∗
ranging from below the critical chirp (−0.35ω∗
c ).
The sub-graph shows the dependence of minimum width on initial re-
duced chirp. For this particular α0, the critical value ω∗
c (red dot),
divides the focusing into bounce-back and crossover regimes. . . . . . . 21

Figure 2.2: Dependence of critical reduced chirp ω∗

a horizontal asymptote of ω∗

c (α0 → ∞) = 2.

c on initial aspect ratio α0, with

. . . . . . . . . . . . . . . 24

Figure 2.3:

t(λ) and integration region in blue. Notice that for λ > λe, we have
t = 0 as Pe is an interior point of those shells

. . . . . . . . . . . . . . 30

Figure 3.1: Schematic representation for the three bunch parameters in SSA . . . . 37

Figure 3.2: Longitudinal width evolution of prolate spheroids with (α0 = 10/75)
focused by diﬀerent initial chirps: (a) 0.7ω∗
c . In
each ﬁgure, red solid line represents the prediction from MFT, dotted
lines represent SSA with diﬀerent emittance ranging from 0 to 0.01µm
or mm · mrads. The SSA with zero emittance and the MFT are in
perfect agreement. Notice that (1) emittance move the waist larger
and earlier comparing SSA and MFT, (2) comparing diﬀerent chirps
in SSA, the bunch evolution driven by critical chirp(b) shows higher
sensitivity for emittance. . . . . . . . . . . . . . . . . . . . . . . . . . . 42

c , (c) 1.5ω∗

c , (b)1.0ω∗

Figure 3.3: Kinetic energy evolution for cross-over case ( −1.5ω∗

c in the longitudinal
direction, corresponding to panel (c) in Fig. 3.2), with solid lines for
MFT and dotted lines for SSA with diﬀerent value of emittance (circle
for KEz and triangle for KEr). The sudden change of direction for
MFT in longitudinal kinetic energy comes from the sign ﬂip of the chirp
as discussed in Sec. 2.

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

x

Figure 4.1: The z − pz phase-space distribution of a warm initial condition. The
initial local momentum ﬂuctuation has a Gaussian proﬁle and the cor-
responding linear chirp is added based on the spatial coordinates of
each electron.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Figure 4.2: Comparison of longitudinal width evolution with diﬀerent initial chirps
from SSA and MD simulations. For each initial chirp, the connected-
dots correspond to the average of the 90 MD simulations, the thick solid
line is the SSA model prediction with the average emittance of the warm
initial conditions. The thick dotted line is the modiﬁed SSA considering
the emittance evolution in Fig. 4.4. The discrepancy between SSA and
MD simulations is driven by the longitudinal emittance decrease, seen
in the simulations, which is conﬁrmed by the good agreement between
the modiﬁed SSA and simulations.

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

Figure 4.3: Comparison of kinetic energy associated with the longitudinal motion
of the electrons relative to the bunch center of mass with diﬀerent
initial chirps from SSA and MD simulations. For each initial chirp, we
have the connected-dots are the corresponding 90 MD simulations, the
thick solid line for SSA prediction with the average emittance of those
warm initial conditions and the thick dotted line for the modiﬁed SSA
considering the emittance evolution in Fig. 4.4.

. . . . . . . . . . . . . 75

Figure 4.4:

(rms) Emittance evolution in (a) transverse εx and (b) longitudinal εz
direction for three diﬀerent initial chirps. . . . . . . . . . . . . . . . . . 76

Figure 4.5: Longitudinal density proﬁle evolution of crossover case (−1.5ω∗

c ) with
the longitudinal width scaled by its standard deviation. In each panel,
solid line represents the average of simulation samples and dotted line
is the density proﬁle for an equivalent uniform spheroid of the same
longitudinal width as the simulation data. The deviation from unifor-
mity is most prominent at crossover. However, tails of the distribution
are present in simulations for all initial chirps. . . . . . . . . . . . . . . 79

Figure 4.6: The variance of the local momentum ﬂuctuations for the case in bounce-

back regime (0.7ωc) on a Log-Log scale.

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

Figure 4.7: The heat transfer from potential to ηi (a) and (b) for bounce-back case
(0.7ωc); (c) and (d) for cross-over case (1.5ωc). The red dotted lines
are the time when the transfer in longitudinal direction switches sign. . 81

xi

Figure 4.8: Longitudinal density proﬁle at tc of one crossover case (−1.5ω∗

c ). The
grey histogram represents the electron distribution along the longitu-
dinal direction of the simulation sample at tc. The green solid line is
the corresponding kernel density estimation (KDE) of the density dis-
tribution. The blue solid line is the corresponding density proﬁle from
the SSA calculation of that sample. The orange dotted line is repre-
senting the density distribution of a uniform spheroid with the same
longitudinal spread. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

Figure 4.9: The longitudinal phase-space (ﬁrst row) and real-space (second row)
distribution evolution of a uniform bunch, which starts with a cylin-
drical shape, zero emittance and a linear longitudinal chirp. Electrons
are colored based on their initial longitudinal position.

. . . . . . . . . 87

Figure 4.10: The variance of the local momentum ﬂuctuations for the cases in the

cross-over regime (1.5ωc) on a Log-Log scale. . . . . . . . . . . . . . . . 89

xii

KEY TO ABBREVIATIONS

Here are some abbreviations and acronyms that appear in this dissertation. The full name

is shown in the ﬁrst appearance of each abbreviation. For a better reading experience

and ease of reference, the list of abbreviations used in this work are summarized below in

alphabetical order.

AG . . . . . . . . Analytic Gaussian

COM. . . . . . Center of mass

DIH . . . . . . . Disorder-induced heating

EOM . . . . . . Equation(s) of motion

FMM. . . . . . Fast multipole method

MD . . . . . . . Molecular dynamics

MFT . . . . . . Mean-ﬁeld theory

KDE . . . . . . Kernel density estimation

K-V . . . . . . . Kapchinsky-Vladimirsky

RF . . . . . . . . Radio-frequency

rms . . . . . . . Root mean square

SSA . . . . . . . Self-similar analytical

UED . . . . . . Ultrafast electron diﬀraction

UEM . . . . . . Ultrafast electron microscopy

xiii

Chapter 1

Introduction

The time scale of fascinating phenomena in the physics world ranges from millions of years

in cosmology to inﬁnitesimally small fractions of seconds in ultrafast phase-transitions

and nuclear reactions. [1, 2] To truly understand the fundamental physics behind these

intriguing phenomena, we need the ability to record the evolution of these systems of

interest at high enough resolution in both the time and space domains. In particular,

resolving the atomic motions throughout phase-transitions leads to the understanding

of critical phenomena in condensed matter physics [3–5]. To fully describe the phase-

transitions with a sub-picosecond time-scale, we need to have the ability to resolve the

temporal resolution at tens of femtoseconds.

1.1 The Ultrafast Electron Microscopy

In practice, many ultrafast phase-transitions are explored experimentally using ultrafast

microscopies with the so-called “pump-probe” scheme [6–11]. As the name suggests,

these experiments consist of two kinds of pulses interacting with the system of interest,

the pump pulse and the probe pulse. A pump pulse comes ﬁrst as the perturbation to the

system of interest. With the excitation from the pump pulse, the system goes through

the phase-transition that we would like to explore. A probe pulse then interacts with

1

the system at a later time (known delay from the pump pulse) to produce a snapshot of

the system at that speciﬁc delay time. This pump and probe process is then repeated

numerous times with various delay times to obtain the time evolution of the system

throughout the phase-transition. In order to resolve these ultrafast phenomena, both the

pump and probe pulses need to be suﬃciently short in time comparing to the time-scale

of the phase-transition process. In addition, for regular ultrafast microscopies, the energy

deposition of the pump and probe pulses are limited so that the system is capable of

going through the same phase-transition and recovering numerous times at a reasonable

repetition rate to allow unfolding the temporal response of the sample.

We are working with the ultrafast research group at Michigan State University, using a

femtosecond laser as the pump pulse to excite the system and the electron bunches as the

probe pulse to measure the system response throughout the phase-transition [12–17]. For

ultrafast electron microscopy (UEM), enough electrons are needed to generate a pattern

with good signal-noise ratio, which translates to roughly 105 to 107 electrons in one single

probe pulse for diﬀraction and 107 to 109 electrons for imaging. [1]

Compressing a bunch with such a high number of electrons introduces technological

hurdles as the space charge eﬀects1 within the probing electron bunch are signiﬁcant at

speciﬁc locations within the microscopy column due to the bunch’s high electron density

at those locations [13–19]. Magnetic lenses are used to mitigate the space-charge eﬀect

induced expansion in transverse direction, while the longitudinal focusing of the electron

bunch is achieved using radio frequency (RF) cavities to mitigate longitudinal space-

charge eﬀects. The schematic design of the beamline is illustrated in Fig. 1.1(a) and

1Here the space charge eﬀect is referring to the strong Coulomb interaction between electrons within

the bunch.

2

Figure 1.1: Ultrafast electron beam column with optics for controlling electron bunch
phase space evolution. (a) Schematic drawing of a prototype UEM beamline. (b)–(d) The
phase space evolution in the injector portion of the beamline (enclosed in dashed line). (e)
The conceptual outline of the atomic grating approach to characterize the energy spread
of the electron bunches. Adapted from “Active control of bright electron beams with RF
optics for femtosecond microscopy” by Williams, J. et. al., 2017, Structural Dynamics,
4(4), 044035.

3

further details can be found in the publications from the Ruan Group [3, 15, 18, 20].

1.2

Importance of our analytical modeling

As the title of this dissertation suggests, we are interested in the analytical modeling

for the evolution of the probing electron bunch, which is critical to the design of next

generation technologies. Simple analytic models are particularly helpful for improving

instrument design.

We focus our analytical study on the uniform ellipsoidal electron bunch, which is

an ideal theoretical object in systems governed by gravitational or Coulomb interactions

[21]. The importance of the uniform ellipsoidal bunch comes from its well-behaved linear

self-ﬁeld, which leads to the maintenance of the uniform charge density proﬁle as the

bunch evolves [22]. In accelerator physics, such a uniform distribution is a prerequisite in

employing techniques such as emittance compensation [23] as well as providing the basis

of other theoretical analyses [17]. It has long been proposed that such a uniform ellipsoid

may be generated through proper control of the transverse proﬁle of a short charged-

particle bunch emitted from a source into vacuum [21, 24–26], and experimental results

have shown that an electron cloud emitted from a photocathode and rapidly accelerated

into the highly-relativistic regime can develop into a ﬁnal ellipsoidal proﬁle characteristic

of a uniform charge distribution [27].

Although we have numerical simulations to illustrate the detailed phase-space distribu-

tion of the electron bunch throughout the evolution, analytical study such as the methods

in this dissertation provides valuable insights about the underlying physics behind the in-

teresting behavior of the bunches.

In addition, the new concepts from the analytical

4

modeling are helpful to better understanding emittance evolution and the broadening of

the energy spread for the probing electron bunches, which is often referred as the ener-

getic Boersch eﬀect [18, 28–30], which is crucial to improve the performance of ultrafast

experiments.

1.3 The Existing analytical models

Currently, the second-order statistics and the corresponding derivative statistics (namely,

the rms emittance) are the key measurements for improving the performance of the prob-

ing electron bunch. The spatial variance in the longitudinal direction is related to the

pulse duration of the electron probe, which is essential to the temporal resolution of the

probe measurement. The variance of longitudinal momentum is another key second-order

statistics, determining the spectral resolution of the probing electron bunch. Therefore,

the existing analytical methods are oriented to present their predictions for the bunch

evolution using the trajectories of those second order statistics.

The diﬀerence in calculation perspective divides the existing models into two cate-

gories: the mean-ﬁeld theory (MFT) and the statistical methods. The MFT here assumes

a uniform density proﬁle with zero emittance and calculates the time derivatives of the

bunch statistics based on that assumption, i.e. the uniform density proﬁle is the cor-

nerstone of the bunch evolution in the MFT. On the contrary, the statistical methods

propagate the electron bunch evolution from the time-derivatives of those second-order

bunch statistics. The force related terms in the time-derivatives of the bunch statistics

are calculated based on the assumption/estimation of the transient density proﬁle, which

is not restricted to one speciﬁc type of function. For example, the density proﬁle can be a

5

uniform bunch for the Kapchinsky-Vladimirsky (K-V) envelope equations or a Gaussian

bunch for the Analytic Gaussian (AG) model. Furthermore, the density proﬁle is not

required to retain the same type during the evolution as the expression of the force is one

of the components to calculate the time derivatives. The driving factor in the statistical

methods is the time derivatives of the bunch statistics, not the speciﬁc density proﬁle.

Here is the brief overview of the three existing models that we will discuss in detail later.

The mean-ﬁeld theory has been utilized in the astrophysics and Coulomb explosion lit-

erature where the mean-ﬁeld eﬀects of a uniform ellipsoidal electron bunch yield ordinary

diﬀerential equations for the bunch statistics. Speciﬁcally, Lin et al. developed a model

for gravitational collapse of a spheroidal ellipsoid that could be written as a system of

diﬀerential equations for the longitudinal and transverse width of the bunch [31]. Similar

techniques regarding the repulsive electrostatic force between the electrons were developed

to model Coulomb explosion from rest [22, 32], which is similar to a time reversed version

of the gravitational collapse. Both techniques require a uniform ellipsoid throughout their

evolution [33].

Michalik and Sipe introduced the Analytic Gaussian (AG) model that predicts not

only the spatial width evolution but also the full phase space evolution [34–36], with the

assumption that the bunch retains its Gaussian density proﬁle throughout the evolution.

As we will discuss in detail later in Sec. 3.3.2, the AG model can be considered as an

intermediate approach between the MFT and the statistical methods. Since AG’s Gaus-

sian density proﬁle uses the three second-order statistics as parameters, the formulation

of the AG model appears to focus on the dynamics of the three bunch statistics while the

original derivation of the AG model is more similar to that of the MFT.

6

As an example of the statistical methods, the well-established Kapchinsky-Vladimirsky

(K-V) envelope equations were initially developed to describe the evolution of uniform

continuous beams [37] in the accelerator physics community. Sacherer provided a sim-

ple perspective which showed that the K-V envelope equations could be derived from

basic statistical considerations with applications of the mean-ﬁeld force from a uniform

distribution [38]. The mathematical form of both the MFT and the AG model can also

be derived from the K-V envelope equations with similar considerations of uniform or

Gaussian density proﬁle, respectively [33].

1.4 Overview

The central focus of this dissertation is to gain more insight into the evolution of the

probing electron bunch from studying the analytical models and numerical simulations.

Here is the detailed structure of this dissertation.

In Chapter 2, we overview the mean-ﬁeld theory (MFT) for studying the compression

process of uniform ellipsoidal electron bunches. MFT provides the expression for electro-

static potential inside uniform and Gaussian ellipsoidal bunches, as well as the evolution

prediction for zero-emittance uniform bunches. C.C. Lin [31] and Grech et. al., [22]

ﬁrst proposed the MFT formalism for uniform bunch evolution from rest. We extend this

model by adding a linear chirp. As a result of this modiﬁcation, we propose a method that

enables researchers to study the crossover phenomena of bunch focusing and its eﬀects on

the evolution of ellipsoidal electron bunches.

In Chapter 3, we use statistical methods to model the evolution of an electron bunch,

particularly a self-similar analytical model. Statistical methods predict the evolution

7

of ellipsoidal bunches through the time derivatives of bunch statistics, which overcomes

the limitation of MFT and the AG model assumptions of the bunch density distribution

proﬁle. We propose the self-similar analytical model, which builds on the MFT and the

Analytic Gaussian model. This new model and its expanded model contributes to the

ﬁeld, as it includes the eﬀect of varying emittance on the evolution of electron bunches. We

also discussed the relationship between the SSA model and the K-V envelope equations.

In Chapter 4, we compare the predictions of the SSA model to numerical simulations

for ellipsoidal bunch compression. In particular, we discuss the impact of varying emit-

tance on bunch evolution. Moreover, we developed a simulation code, using it in the

electron bunch compression process.

Chapter 5 presents the conclusions and discussion. In particular, we discuss the direc-

tions for future research, such as developing analytic models for emittance growth, given

the importance of these two concepts on bunch evolution.

8

Chapter 2

Mean-ﬁeld theory for dynamics of

ellipsoidal electron bunches

In this chapter, we propose a mean-ﬁeld theory (MFT) to describe evolution of an el-

lipsoidal probing electron bunch in ultrafast experiments. We ﬁrst review the mean-ﬁeld

theory for uniform ellipsoidal bunches and then we discuss the evolution of a Gaussian

bunch. We conclude this chapter by discussing the underlying assumptions of the MFT

such as the linearity assumption and self-similar evolution.

We designated the name “MFT” here speciﬁcally to the formalism that studies uni-

formly charged ellipsoids with perfect momentum-position correlation (zero emittance).

This formalism takes advantage of the linear ﬁeld calculated using mean-ﬁeld theory and

the linearity of the uniform spheroidal bunch.

2.1 Mean-ﬁeld Theory for a uniform spheroidal bunch

For modeling electron bunch evolution, the majority of existing UEM literature discusses

uniform distributions, which stems from general potential theories regarding both gravi-

tational [31] and Coulomb interactions [21, 22]. The advantage of working with uniform

bunch dynamics is that they evolve self-similarly.

In other words, the uniform bunch

9

stays uniform throughout the evolution. This self-similar evolution greatly simpliﬁes the

modeling process and the bunch evolution can be captured by tracking the bunch size

parameter in each direction. The self-similar evolution comes from the linear self ﬁeld

inside the uniform bunch. As Lin et. al, pointed out in their paper for gravitational

collapse [31], the essence is that the potential within a uniform spheroid is quadratic in

Cartesian coordinates, i.e., the self ﬁeld linearly depends on the spatial position. For

Coulomb interactions, in electron bunches, a similar idea was applied to study Coulomb

explosion from rest with zero emittance and achieved good agreement with numerical

simulations [22, 32].

To provide necessary context for understanding the MFT, we ﬁrst review the non-

relativistic MFT for a uniform spheroidal bunch, and demonstrate how linear initial

chirps, i.e. the linear correlation between momentum and position, can be added into

the model to describe the bunch focusing process. Here, the bunch focusing is referring

to the temporal compression of the probing electron bunch by an radio-frequency (RF)

cavity, or a transverse compression by a magnetic focusing lens along the beamline. Also,

details about the initial conditions are examined, and the discussion naturally leads to

the identiﬁcation of a critical chirp that demarcates two qualitatively diﬀerent regimes of

bunch behavior according to MFT: “bounce-back” where electrons reverse their motion

and never cross the bunch center and “cross-over” where electrons cross the bunch center.

2.1.1 Mean-ﬁeld theory for a uniform ellipsoidal bunch

The mean-ﬁeld theory we are referring to here, treats the electron bunch as an ideal

homogeneous continuum object. This simpliﬁcation can be considered as taking the con-

10

tinuum limit of the actual bunch consisting of a ﬁnite number of discrete electrons in the

experiments and simulations, i.e. MFT is essentially modeling a bunch with N identical
particles where N → ∞ and each particle carries charge amount Q/N , with Q being the

total charge of the entire bunch. Thanks to this simpliﬁcation, all the nice features of

mean-ﬁeld theory are applicable to the uniform electron bunch, while the limitation of this

simpliﬁcation will be discussed later. The following derivation of the Coulomb potential

and ﬁeld for uniform ellipsoids are generally based on Chapter 2 in MacMillan’s book [39]

for the gravitational interaction.

The surface of a given ellipsoid with semi-axes (a, b, c) is deﬁned by equation:

ξ2
a2 +

η2
b2 +

ζ2
c2 = 1

(2.1)

Let the interior point for which the potential is to be computed be P (x, y, z). On taking

P as the origin of a spherical coordinates system (ρ, ϕ, θ) with the transformation:

ξ = x + ρ cos ϕ cos θ

η = y + ρ cos ϕ sin θ

ζ = z + ρ sin ϕ

(2.2a)

(2.2b)

(2.2c)

and the corresponding charge element

dq = ρcρ2 cos ϕdϕdθdρ

(2.3)

with the volume charge density ρc = n · e (typically n is the electron number density).

11

Then, the electrostatic potential at point P (with respect to the potential zero point at

inﬁnity) is expressed as:

V (x, y, z) =

(cid:90)

E

ρc
4πε0

dq
ρ

= κ

(cid:90) + π

2

− π
2

(cid:90) 2π

(cid:90) ρ1(θ,ϕ)

0

0

We substitute κ = ρc
4πε0

to simplify the derivation.

ρ cos ϕdϕdθdρ

(2.4)

The upper limit ρ1(θ, ϕ) of the integration with respect to ρ is a function of θ and ϕ,

since the integration is from P to a point on the surface of the ellipsoid. So, we can insert

Eq. 2.2a into Eq. 2.1 for ρ1 as:

Aρ2

1 + 2Bρ1 + C = 0

where

A =

B =

C =

+

+

cos2 ϕ sin2 θ

b2

+

sin2 ϕ

c2

y cos ϕ sin θ

b2

+

z sin ϕ

c2

cos2 ϕ cos2 θ

a2

x cos ϕ cos θ

a2
y2
b2 +

x2
a2 +

z2
c2 − 1

After a lengthy derivation,1 we can have

(cid:18)cos2 ϕ cos2 θ
(cid:90) + π
(cid:90) 2π
(cid:90) 2π
(cid:90) + π

− π
2

cos ϕdϕdθ

0

2

a2

− π
2

0

A

2

V = κ

− κ
2

C

· x2
a2 +

cos2 ϕ sin2 θ

b2

· y2
b2 +

sin2 ϕ

c2

· z2
c2

(2.5)

(2.6a)

(2.6b)

(2.6c)

(cid:19) cos ϕdϕdθ

A2

(2.7)

1The detailed derivation can be found in the Appendix

12

To simplify the above expression, one may use the following:

(cid:90) + π

2

(cid:90) 2π

− π
2

0

W =

κ
2

cos ϕdϕdθ

A

(2.8)

The simpliﬁed form of the potential of an interior point then reduces to:

V =

=

1
a

(cid:18) 1

∂W
∂a
∂W
∂a

a

x2 +

1
b
− W
a2

∂W
∂b

(cid:19)

y2 +

x2 +

1
c

(cid:18)1

b

∂W
∂c
∂W
∂b

(cid:19)

z2 − CW
− W
b2

y2 +

(cid:18)1

c

∂W
∂c

− W
c2

(cid:19)

(2.9)

z2 + W

Since W is a function of the semi-axes (a, b, c), then so are all of its derivatives. Therefore,

the coeﬃcients of the quadratic term x2, y2 and z2 are functions of a, b, and c only. Also,

W is the potential at the center of the ellipsoid if taking P (x, y, z) = (0, 0, 0). W can be

further reduced to the following form,2

W = πκabc

.

(2.10)

Thus the term related to the derivative of W with respect to a reads:

(cid:90) ∞

0

ds

(cid:112)(a2 + s)(b2 + s)(c2 + s)

(cid:112)(a2 + s)(b2 + s)(c2 + s)

1√
a2 + s

− 1
2

−

ds

∂a(a2 + s)
(cid:90) ∞
(cid:90) ∞

∂a

0

1

a2 + s

0

1
a

∂W
∂a

− W

a2 =

πκabc

a2

= −πκabc

(cid:112)(b2 + s)(c2 + s)

ds

(2.11)

Therefore, the potential inside of a uniform ellipsoidal electron bunch may be expressed

2The detailed derivation can be found in the Appendix

13

as:

V (x, y, z) =

ρc
4ε0

abc

(cid:90) ∞

0

(cid:18)
1 − x2
a2 + s

− y2
b2 + s

− z2
c2 + s

(cid:19)

(cid:112)(a2 + s)(b2 + s)(c2 + s)

ds

(2.12)

with the zero point of potential at inﬁnity.

The corresponding electrostatic ﬁeld at the interior point P (x, y, z) is,

(cid:126)E(x, y, z) = Ex ˆx + Ey ˆy + Ez ˆz

(2.13)

with ˆx, ˆy, and ˆz representing the unit vectors, respectively. The components of the ﬁeld

can be represented by the following equations:

Ex(x, y, z) = x · ρcabc
2ε0
Ey(x, y, z) = y · ρcabc
2ε0
Ez(x, y, z) = z · ρcabc
2ε0

(cid:90) ∞
(cid:90) ∞
(cid:90) ∞

0

0

0

ds

(a2 + s)(cid:112)(a2 + s)(b2 + s)(c2 + s)
(b2 + s)(cid:112)(a2 + s)(b2 + s)(c2 + s)
(c2 + s)(cid:112)(a2 + s)(b2 + s)(c2 + s)

ds

ds

(2.14a)

(2.14b)

(2.14c)

Speciﬁcally for spheroids, we introduce the radial coordinate, r, to take advantage of

the rotational symmetry of spheroids. The radial coordinate r is expressed as:

(cid:113)

r =

x2 + y2

(2.15)

Although the detailed calculations below are for prolate spheroids (a = b < c), similar

results are valid for general, uniformly charged ellipsoidal bunches with arbitrary semi-

axes (a, b, c). We can then rewrite the electrostatic ﬁeld (Eq. 2.14) as a linear function of

14

the cylindrical coordinates, as follows:

(cid:126)E(r, z) = Er(r)ˆr + Ez(z)ˆz

(2.16)

with ˆr and ˆz representing the radial and longitudinal unit vectors, respectively, and

Er(r) =

Ez(z) =

ξr(α) · r

ξz(α) · z

ρc
2ε0
ρc
2ε0

where the spheroid aspect ratio α is deﬁned as:

α =

a
c

and where the corresponding geometry coeﬃcients are

(cid:90) ∞
(cid:90) ∞

0

ξr(α) = α2

ξz(α) = α2

ds

(α2 + s)2(1 + s)1/2

ds

0

(α2 + s)(1 + s)3/2

(2.17a)

(2.17b)

(2.18)

(2.19a)

(2.19b)

Based on Eq. 2.17 to Eq. 2.19, we can get the equations of motion through the aspect

ratio of the bunch, and the entire bunch evolution history can be obtained with proper

initial conditions.

15

2.1.2 Self-Similar Evolution and Equations of Motion

Before we discuss the equations of motion for the bunch evolution, we need to discuss the

self-similar assumption that enables this simple MFT formalism. The linear dependence of

the electrostatic ﬁeld on spatial position results in the preservation of bunch uniformity

provided that the initial momentum-position proﬁle does not start out non-linear.

In

other words, the charge density must be homogeneous across the bunch, while the bunch

size and aspect ratio evolves. This self-similar evolution greatly simpliﬁes our analysis as

the formulation presented in preceding sections always applies to the bunch throughout

the process, and the evolution is reduced to the determination of two degrees of freedom.

Speciﬁcally, the temporal evolution of the entire bunch can be represented by the evolution

of two unit-less scaling functions, R(t) and Z(t). Since these formulas hold, the trajectory

of any particle with initial positions (r0, z0) inside the uniform spheroid at a later time t

can be expressed as by (r0R(t), z0Z(t)), where R and Z are independent of the arbitrary

initial position (r0, z0). Thus, the parameters for describing the bunch change according

to:

• the semi-axes of the spheroids as (a, c) = (a0R, c0Z),

• the transient aspect ratio as α(t) = α0 · R
Z ,

• the charge density from conservation of charge3 as ρc(t) =

ρc0
R2Z

=

n0e
R2Z

,

• the longitudinal width as σ2

z (t) = σ2

z0 · Z2, with σ2

z = c2
5 .

In the non-relativistic limit, each individual electron [j] with initial position (r0, z0)
3Ntotal = ρc0 · (4π/3)a2

0c0 = ρc(t) · (4π/3)a2c

16

follows the Newton’s second law as:

d2
dt2 (r0R) =
d2
dt2 (z0Z) =

eEr(α)

me

eEz(α)

me

=

=

n0e2
2ε0me
n0e2
2ε0me

ξr(α)
R2Z
ξz(α)
R2Z

· (r0R)
· (z0Z)

(2.20a)

(2.20b)

Similarly, the bunch evolution then follows the equations of motion (EOM) in such an

electric ﬁeld, which can be reduced to the following two dimensionless ordinary diﬀerential

equations (ODEs):

with unit-less reduced time τ ,

d2R
dτ 2 =
d2Z
dτ 2 =

ξr(α)
RZ
ξz(α)
R2

(cid:115)

τ = t ·

e2n0
2ε0me

= t · Ω0

(2.21a)

(2.21b)

(2.22)

initial electron number density n0, and electron mass me. Notice that ﬁrst, the time

(cid:114)

scaling factor Ω0 =

ωp0√
2

, where ωp0(n0) =

e2n0
ε0me

is the initial plasma frequency. Second,

bear in mind that the geometry coeﬃcients ξr and ξz are solely dependent on the transient

aspect ratio α rather than speciﬁc values of a and c.

This formula shows that starting with the same initial conditions for the ODE, bunches

with the same initial aspect ratio α0 but with diﬀerent initial densities n0, will have

behaviors which only diﬀer by the time scaling factor Ω0 determined by the initial density

n0. The ones with higher density evolve faster, but the evolution trajectory reduces to

17

the one with lower density after the evolution time is scaled according to the initial Ω0.

2.1.3

Initial condition and Coulomb Explosion

We present the initial condition for solving uniform electron bunch evolution with the

ODE described above. In order to achieve this, we rely on one particular example from

Grech et al., [22], which begins with the initial conditions corresponding to a uniform

spheroidal electron bunch expanding under the Coulomb force from rest, a phenomenon

often referred to as Coulomb explosion.

According to the deﬁnition of the two scaling functions, for the Coulomb explosion,

we have:

R(τ = 0) = 1

Z(τ = 0) = 1

(cid:12)(cid:12)(cid:12)(cid:12)τ =0
(cid:12)(cid:12)(cid:12)(cid:12)τ =0

dR
dτ

dZ
dτ

= 0

= 0

(2.23a)

(2.23b)

(2.23c)

(2.23d)

Speciﬁcally, Eq. 2.23a and Eq. 2.23b represent the initial scaling of the spheroid and are

by deﬁnition equal to 1. Meanwhile, Eq. 2.23c and Eq. 2.23d represent the initial change

rate of the scaling function R and Z, or the initial velocity of the expansion. As the bunch

starts from rest in this Coulomb explosion calculation, these rates equal to zero.

Using these initial conditions, the predictions derived from the ODEs (Eq. 2.21) are

found [22, 32] to be in good agreement with molecular dynamics (MD) simulations for

time-dependent energy distribution and particle-in-cell (PIC) simulations for temporal

18

spheroid radii evolution.

2.1.4 Critical Chirp ω∗

c , Crossover and Bounce-back

In addition to setting the initial velocity of the scaling faction (Eq. 2.23c and Eq. 2.23d)

to zero for bunches starting from rest, we can also assign non-zero initial velocities,4 for

bunches start with initial linear chirp. For example, we can have initial velocities as:

(cid:12)(cid:12)(cid:12)(cid:12)τ =0
(cid:12)(cid:12)(cid:12)(cid:12)τ =0

dR
dτ

dZ
dτ

= 0

= ω∗

(2.24a)

(2.24b)

where ω∗ is proportional to a longitudinal linear chirp. By linear chirp, we mean a

linear correlation between momentum and position. Speciﬁcally, the linear chirp in MFT

means that we are assuming for each individual particle, the longitudinal momentum pz

is proportional to its longitudinal position z.

The addition of this longitudinal linear chirp enables MFT to model bunches with

a linear chirp in the longitudinal direction, such as the probing electron bunch after a

longitudinal focusing lens, e.g. an RF cavity. As such, the initial momentum in the

longitudinal direction of any electron [j] with its longitudinal position z[j] is expressed

as:

pz[j] = Cz · z[j] = me · (ω∗Ω0 · z[j])

(2.25)

where Cz is the actual linear chirp in simulations and experiments. We refer to ω∗ as

reduced linear chirp because it serves as a linear chirp in the dimensionless EOM and

4Instead of being any function of r or z

19

scales from the actual chirp in experiments only by Ω0 as a result of the diﬀerence in

the reduced time τ . If the initial reduced chirp is negative, Z will initially decrease, i.e.

the bunch will be compressed in the longitudinal direction in the initial evolution. As
ω∗ is proportional to the actual chirp via a parameter entirely dependent on density, we
drop the “reduced” from ω∗ when we refer to only the reduced chirp, but we keep the

“reduced” when we are comparing the reduced chirp in the model with the actual chirp

from the experiments and simulations.

The initial bunch density determines the time scale of the evolution via Ω0. Therefore,

the categorization of the focusing process of a uniformly charged spheroid is entirely

determined by the initial reduced chirp and the initial aspect ratio. The eﬀect of aspect

ratio on bunch evolution has been well studied [22]. Therefore, we examine the impact

of the reduced chirp on the bunch focusing process at a pre-speciﬁed initial aspect ratio,

α0 = 10

75, which is typical for .

We deﬁne the critical time, τc, as the time when the bunch reaches its minimum width,

min(σ2

z ) during the focusing process. As Eq. 2.21b is only dependent on 1

R2 , which is
ﬁnite throughout our investigation, Z can become zero at some point in the compression

process if the initial chirp provides suﬃcient focusing power to overcome the repulsive

Coulomb force, for zero emittance bunch. We call the smallest magnitude of the initial

reduced chirp that is able to compress Z to zero (i.e. min(σ2
chirp, ω∗

c . This critical reduced chirp is also unitless, just as all the other quantities in

z ) = 0) the critical reduced

the EOM and the initial condition. That means uniform spheroidal bunches with the

same aspect ratio share the same reduced critical chirp as we will discuss later in this

section. The focusing process then falls into two categories characterized by whether

20

min(σ2

z ) reaches zero. As shown in Fig. 2.1, for ω∗ < ω∗

c , the longitudinal width of the

Figure 2.1: Longitudinal width evolution Z(τ ) of prolate spheroids with (α0 = 10/75)
driven by diﬀerent initial chirps in numeric solutions of MFT ranging from below the
critical chirp (−0.35ω∗
c ). The sub-graph shows the dependence of
minimum width on initial reduced chirp. For this particular α0, the critical value ω∗
c (red
dot), divides the focusing into bounce-back and crossover regimes.

c ) to well above (−2.0ω∗

bunch is never zero. However, there is a time at τc, at which the bunch size reaches a
minimum during the compression process. When τc = 0 then ω∗ = 0, and it increases
with ω∗ increasing. This trend continues until the critical reduced chirp is reached where

the minimum width reaches zero and the corresponding critical time τc goes to inﬁnity
(τc → ∞). When compressed by a chirp that is larger than the critical chirp, the bunch

will overcome the Coulomb repulsion and be compressed through a longitudinal crossover,

21

024681012140.00.20.40.60.81.0Z0.35 c0.70 c1.00 c1.50 c2.00 c0.00.20.40.60.00.20.40.60.81.0Min(Z)cbounce-backcrossoveras electrons starting from one side of the bunch cross the center of mass and then begin

to continue moving on the other side without changing direction of motion. We refer

to this type of process as the “cross-over” regime. In this regime, further increasing the

chirp will decrease the critical time for the bunch to reach its minimum width, but the

minimum width of the bunch during such compression will always be zero. In contrast,

we call the regime below critical chirp the “bounce-back” regime as the particles follow

trajectories that start out towards the center of mass and reverse their direction after the

bunch reaches the minimum width, without any particle crossing.

The linearity of MFT (linear initial condition and linear self-ﬁeld) indicates that for

a bunch in the cross-over regime, all crossover incidents happen simultaneously across

the entire bunch at the critical time (τc), creating a singularity in the EOM where the

longitudinal width goes to zero. To work through this singularity while numerically solving

the EOM, we use a small time step to propagate the EOM until Z falls below zero. Then,

we stop the calculation and ﬂip the value of both the longitudinal position scaling, Z, and

the longitudinal momentum, pz to positive. The same EOM are then used to integrate

the parameters for the expansion after crossover. In eﬀect, this two-step process skips

the singularity by inﬁnitesimal step size in time. This two-step process is equivalent to

assuming that the change of the bunch in the transverse direction is negligible during

the inﬁnitesimal time step around τc, which is validated by the fact that the transverse

change rate or velocity is not divergent throughout the crossover. In addition, this two-

step process also implies that in the cases where Z passes through zero, momentum does

not change sign.

Analogous to the linear chirp we added in longitudinal direction, a radial chirp can also

22

be introduced by altering the radial initial velocity in Eq. 2.23c to a non-zero value. The

linear chirp in the radial direction may be combined with the longitudinal chirp to model

bunches compressed in both degrees of freedom simultaneously. However, in contrast to

the longitudinal direction, MFT predicts that there is no such critical chirp in the radial

direction. This occurs because the change in radial velocity is proportional to the inverse

of R from Eq. 2.21, which diverges as the bunch is compressed radially, preventing R

from reaching zero. Compressing the bunch radially is equivalent to compressing the

ˆx and ˆy simultaneously, and in this case there is only the bounce-back regime in the

radial direction. Here, we focus on the longitudinal focusing where both bounce-back and

cross-over regimes are accessible in MFT.

The (longitudinal) critical reduced chirp is an interesting concept that is important

but new to the community.

It is important because the bunch evolution is sensitive

to emittance, especially when compressed by the critical chirp as we will show in next
chapter. One important feature of the critical reduced chirp ω∗

c , is its exclusive dependence

on initial aspect ratio α0, which stems from the governing EOM solely depending on the
aspect ratio. In Fig. 2.2, we present the reduced critical chirp ω∗
c as a function of the
initial aspect ratio α0. Speciﬁcally, note that for large α0, we will always have α (cid:29) 1,

where the geometry coeﬃcients ξr and ξz in the EOM can be approximated in closed

forms as [22]:

ξr(α → ∞) (cid:39) π
2α
ξz(α → ∞) (cid:39) 2 −

→ 0

πα2

(α2 − 1)3/2

→ 2

(2.26a)

(2.26b)

23

Figure 2.2: Dependence of critical reduced chirp ω∗
horizontal asymptote of ω∗

c (α0 → ∞) = 2.

c on initial aspect ratio α0, with a

The EOM reduces to the following:

R (cid:39) 1

d2Z
dτ 2 = 2

(2.27a)

(2.27b)

In other words, the trajectory of Z reduces to a constant acceleration. Therefore, the

corresponding critical reduced chirp for large initial aspect ratio is simply:

Z(t) = Z0 − ω∗

c t +

ξz(∞)t2

1
2

⇒

ω∗
c = 2

(2.28)

24

10210110010110210310400.010.101.002.00c0=10/75as we can see from Fig. 2.2.

2.2 Coulomb Potential inside a Gaussian Bunch

In addition to uniform bunches, researchers have been studying Gaussian density proﬁles

of electron bunches [21] for various reasons. Gaussian bunches are appealing because the

laser pulses generating the electrons from photo-cathodes generally have Gaussian pro-

ﬁles. In addition, the rich techniques available for Gaussian integrals provide a promising

path toward a closed-form analytical model, which is also appealing to the theoretical

community.

For the Gaussian bunch we are interested in, we assume that the charge density at

point (x, y, z) depends on its scaling variable λ [40]:

λ(x, y, z) =

x2
a2 +

y2
b2 +

z2
c2

(2.29)

where the semi-axes (a, b, c) equal the standard deviation of the entire bunch in each

direction (σx, σy, σz). Then, the charge density proﬁle can be expressed as a function of

λ:

ρ(λ) = ρ(

x2
a2 +

y2
b2 +

z2
c2 ) =

(cid:90)

dλ · f (λ) · δ(λ − x2

a2 − y2

b2 − z2
c2 )

Accordingly, a Gaussian bunch can be modeled with:

fG(λ) =

Q

π3/2abc

e−λ

where Q is the total charge of the bunch. In other words, this bunch is constructed by a

25

(2.30)

(2.31)

series of so-called homogeneous “scaled” shells, where the shells are bounded by ellipsoids

with varying λ and the local volume charge density is the same within each shell with the

value fG(λ).

We start with the contribution to the Coulomb potential from one single shell. We

then integrate the contributions of each shell over the entire bunch yielding the Coulomb

potential and ﬁeld at any interior point of a Gaussian bunch. To conclude, we discuss the

numerical strategy to model the Gaussian bunch evolution.

2.2.1 Potential Contribution from one single shell

We already have the potential for interior points P (x, y, z) of a uniform ellipsoid of semi-

axes (a, b, c) with density ρc from Eq. 2.12. Now, consider a similar ellipsoid of the

√
same volume charge density ρc with scaled semi-axes (
λa,

√

√

λb,

λc). Accordingly, the

surface of this scaled ellipsoid is deﬁned by:

x2
a2 +

y2
b2 +

z2
c2 = λ

We refer to this kind of arrangement as “similar and similarly placed” ellipsoids. Then,

if Ex2 is the x-component of the self ﬁeld at P (x, y, z) then,

Ex2(x, y, z) = x · ρcabc
2ε0

(cid:90) ∞

0

(λa2 + s)(cid:112)(λa2 + s)(λb2 + s)(λc2 + s)

λ3/2ds

.

(2.32)

26

In this case, the integration variable can be substituted with τ = s/λ, yielding the same

electrostatic ﬁeld as in Eq. 2.17:

Ex2(x, y, z) = x · ρcabc
2ε0

(a2 + τ )(cid:112)(a2 + τ )(b2 + τ )(c2 + τ )

dτ

= Ex(x, y, z)

(2.33)

(cid:90) ∞

0

and similarly, for Ey2 and Ez2. Therefore, the ﬁeld at any point of the vacuum space inside

such homogeneous shell, which is bounded by two scaled and similarly placed ellipsoids,

is zero. Meanwhile, if λ > 1, the potential in the hollow space has the constant value
V = (λ − 1)W , where W is deﬁned in Eq. 2.8.

If λ is very close but just slightly larger than 1, then the shell can be considered as a

thin shell. Although the charge volume density is homogeneous within this thin shell, the

surface charge density is not uniform, because the thickness is not a constant across the

thin shell. The thickness of any point of the shell is proportional to the distance h from

the center to the tangent plane at that point [41], which can be expressed as:

h(x, y, z) =

(cid:114)
a2 + y2
x2
a4 + y2
x2

b2 + z2
c2
b4 + z2
c4

(cid:114)

=

1

a4 + y2
x2

b4 + z2
c4

Therefore, the surface charge density on the shell reads:

ρc ·(cid:104)

σ(x, y, z) = lim
λ→1

(λ1/2 − 1) · h(x, y, z)

(cid:105)

(cid:114)

= lim
dλ→0

1

2 · ρc · dλ
a4 + y2
x2

b4 + z2
c4

(2.34)

(2.35)

Meanwhile, the surface charge density for a charged ellipsoidal conductor with the same

27

semi-axes (a, b, c) shares the same form as:

(cid:18) x2

(cid:19)− 1

2

y2
b4 +

z2
c4

(2.36)

σ =

Q

4πabc

a4 +

with Q being the total charge on the conductor [41, 42]. The similarity in the surface

charge distribution suggests that5 we can use the existing potential and ﬁeld expressions

of a charged ellipsoidal conductor for our homogeneously charged thin shell. With a total

charge Qλ on the homogeneous thin shell λ in our Gaussian bunch6 reads:

Qλ = lim
dλ→0

ρ(λ)dV = 2πabc ·

ρ(λ)λ1/2dλ

lim
dλ→0

(2.37)

Therefore, the contribution from any thin shell λ to the Coulomb potential at a point

of interest Pe(xe, ye, ze) can be modeled as:



(cid:82) ∞
(cid:82) ∞

ξ

0

dξ

Rλ(ξ)

dξ

Rλ(ξ)

λe ≥ λ

λe < λ

(2.38)

1

4πε0

1

4πε0

Qλ
2

Qλ
2

φλ(xe, ye, ze) =

with the corresponding Rλ(ξ)

(cid:113)

Rλ(ξ) =

(λa2 + ξ)(λb2 + ξ)(λc2 + ξ) = λ3/2R(

ξ
λ

)

(2.39)

and the lower limit of integration in ξ for exterior point cases comes from the largest root

5In fact, the surface charge density on the ellipsoidal conductor is proportional to the fourth root of

the total curvature of the surface. A detailed discussion can be found in [41]

6We are using the scaling variable λ of the inner surface of the shell from Eq. 2.29 to designate each

thin shell across our Gaussian bunch

28

of the following equation:

x2
e

λa2 + ξ

+

y2
e

λb2 + ξ

+

z2
e

λc2 + ξ

= 1

In other words, when the scaling variable λe of Pe:

λe = λ(xe, ye, ze) =

x2
e
a2 +

y2
e
b2 +

z2
e
c2

(2.40)

(2.41)

is larger than λ, the scaling factor of the thin shell, Pe is an exterior point and the lower

limit of integration is the corresponding ξ from Eq. 2.40. While if λe is smaller than λ,

then Pe is an interior point and the lower limit of integration is zero.

2.2.2 Potential at an interior point of a Gaussian bunch

The potential at Pe due to the entire Gaussian bunch can be obtained by integrating the

corresponding contribution over all the thin shells, as:

φ(xe, ye, ze) =

φλ<λe +

φλ>λe

(cid:88)

=

=

=

abc
4ε0
abc
4ε0
abc
4ε0

(cid:88)

(cid:90) λe
(cid:90) ∞
(cid:90) ∞

0

0

ρ(λ)λ1/2dλ

ρ(λ)λ1/2dλ

(cid:90) ∞

ρ(λ)dλ

(cid:90) ∞
(cid:90) ∞

ξ(λ)

(cid:90) ∞

λe

+

abc
4ε0

dξ

Rλ(ξ)
dξ

(cid:90) ∞

dξ

0

Rλ(ξ)

ρ(λ)λ1/2dλ

ξ(λ)

λ3/2R(ξ/λ)

dt

R(t)

(2.42)

0

t(λ)=ξ(λ)/λ

29

Here the relationship between t and λ for point Pe can be expressed as:

λ(xe, ye, ze; t) =

x2
e

a2 + t

+

y2
e

b2 + t

+

z2
e

c2 + t

(2.43)

derived from Eq. 2.40. The corresponding region of integration is shown as the colored

area in Fig. 2.3.

Figure 2.3: t(λ) and integration region in blue. Notice that for λ > λe, we have t = 0 as
Pe is an interior point of those shells

In order to take the advantage of our knowledge about the charge distribution function

ρ(λ), we interchange the order of integration:

(cid:90) ∞

0

φ(xe, ye, ze) =

abc
4ε0

ρ(λ)dλ

(cid:90) ∞

dt

ξ(λ)/λ

R(t)

(cid:90) ∞

0

(cid:90) ∞

λ(t)

dt

R(t)

=

abc
4ε0

ρ(λ)dλ

(2.44)

assuming the potential is zero at inﬁnity. If the potential is normalized such that the zero

30

et=()/t()point of the potential is at the center of the ellipsoid, then the potential can be written

as:

φ(xe, ye, ze) = − abc
4ε0

dt

0

R(t)

0

(cid:90) ∞

(cid:90) λ(t)

ρ(λ)dλ

(2.45)

For our Gaussian bunch, we insert the charge density proﬁle as Eq. 2.31:

ρ(λ) =

Q

π3/2abc

e−λ

(2.46)

into Eq. 2.44. Therefore, the potential expression for an interior point of a Gaussian

bunch reads:

φ(xe, ye, ze) =

=

(cid:90) ∞
(cid:90) ∞

0

exp

(cid:90) ∞

λ(t)

dt

(cid:20)
R(t)
− x2
a2 + t

e

e−λdλ

− y2
b2 + t

e

− z2
c2 + t

e

(cid:21) dt

R(t)

(2.47)

abc
4ε0

Q

π3/2abc
Q

4π3/2ε0

0

This expression is consistent with the existing literature [40]. In addition, the potential

on an interior point of a uniform bunch in Eq. 2.12 can be reproduced as we take ρ(λ) to

a constant ρc in Eq. 2.45, with a trivial diﬀerence in the zero point of Coulomb potential.

Since the only assumption involved is the ellipsoidal symmetry of the density proﬁle (Eq.

2.30) with no restriction on the speciﬁc form of the density proﬁle, this potential expression

(Eq. 2.44) is valid for bunches with same or higher symmetry like spherical (a = b = c)
and spheroidal (a = b (cid:54)= c) symmetry.

Since the resulting electrostatic ﬁeld from the Gaussian bunch is non-linear, it is

not reasonable to apply the MFT formalism presented in the previous section for the

evolution of a Gaussian bunch. A detailed discussion about this linearity assumption

will be presented in the next section. Nevertheless, this potential and ﬁeld formula from

31

MFT will be useful for our study when we apply statistical methods for the evolution of

Gaussian bunches in the following chapter.

2.3 Discussion about the linearity assumption of MFT

The key assumption of the MFT formalism (Eq. 2.21 as EOM with linear initial condition)

is the linearity of the uniform bunch. This fundamental assumption stems from the linear

self-ﬁeld of uniform ellipsoidal bunches. The linearity of the entire system is maintained as

long as the initial condition is also linear. Therefore, for MFT, the linearity assumption is

equivalent to the self-similar evolution assumption, which greatly simpliﬁes the expression

of the uniform bunch evolution.

This linearity assumption is also one of the major limitations to MFT. Because of

this, MFT cannot be applied to non-linear systems, such as the Gaussian bunch. Even

if we divide the Gaussian bunch into scaled thin shells and treat them individually, the

Gaussian bunch still drastically violates the linearity/self-similar evolution assumption of

MFT, because its self-ﬁeld is non-linear. Further complicating the situation, the non-

linear ﬁeld leads to diﬀerent aspect ratios for shells across the bunch, which then breaks

our scaled thin shell set-up for potential calculation in MFT.

In addition, the MFT also assumes zero-emittance of the bunch, which leads to the

singularity at the cross-over point for cases in the cross-over regime. This singularity

predicted by the MFT will disappear with the presence of a non-zero emittance. Speciﬁ-

cally, the non-zero local momentum ﬂuctuation, which cannot be completely eliminated,

leads to a non-zero minimum width at the cross-over point (tc) for the cross-over cases,

comparing to the zero minimum width predicted by the MFT. Since emittance is one of

32

the key factors for bunch performance, we need to study the emittance eﬀect on bunch

evolution carefully with the help of other methods.

Despite these limitations, MFT can be still be used to predict bunch evolution even

when its assumptions are not strictly satisﬁed. When the assumptions are not fully met,

the predictions will deviate from the actual evolution. The extent of this discrepancy

depends on how far the bunch initial condition is from a uniform bunch with a perfect

linear chirp used in MFT. A detailed discussion about the impact of diﬀerent kinds of

assumptions will be presented in the next chapter.

To reduce the limitations associated with linearity and zero-emittance, we turn to

statistical methods, which focus on the evolution of bunch statistics rather than the exact

form of the density proﬁle.

33

Chapter 3

Statistical Methods

In this chapter, we present statistical methods to model the evolution of electron bunches.

We begin with an overview of the self-similar analytical (SSA) model for uniform el-

lipsoidal bunches. We then discuss the concept of emittance and its impact on bunch

evolution by comparing SSA to MFT. To follow, we discuss a semi-statistical model for

Gaussian bunches — the Analytic Gaussian (AG) model — to further illustrate the un-

derlying assumptions of these statistical methods. We conclude this chapter by discussing

the advantages and limitations of AG, SSA, and the equivalent Kapchinsky-Vladimirsky

(KV) envelope equations.

3.1 The Self-similar Analytical Model

The SSA model, and other statistical methods, extract the bunch evolution trajectory

directly from the time-derivatives of corresponding bunch statistics. We refer to these

types of models as statistical methods.

In contrast, the uniform bunch evolution in

the MFT is found by integrating the uniform density proﬁle with corresponding time-

dependent parameters.

To properly introduce SSA and other statistical methods, we ﬁrst discuss the dynamics

of bunch statistics for a uniform ellipsoidal bunch. In each degree of freedom (i = T, z,

34

where T and z represents the transverse and longitudinal direction, respectively), there

are three second order statistics in the position-momentum phase space:1 1) the variance

of the position σ2

i , 2) the variance of the momentum σ2
pi

, and 3) the correlation between

position and momentum σi,pi

. Speciﬁcally, in the center of mass frame, we can have these

three statistics:

i = (cid:104)i2(cid:105)
σ2
= (cid:104)p2
i(cid:105)
= (cid:104)ipi(cid:105)

σ2
pi

(3.1a)

(3.1b)

(3.1c)

where i = x, y,or z and the (cid:104)(cid:105) operator yields the mean, e.g. (cid:104)zpz(cid:105) =(cid:82)(cid:82) zpzf (r, p)drdp

σi,pi

for a continuous distribution in analytical models.

The non-relativistic evolution of the position and momentum can be written as:

dt + O(dt2)
i(t + dt) = i(t) +
pi(t + dt) = pi(t) + Fi(t)dt + O(dt2)

pi(t)
m

(3.2a)

(3.2b)

with the force distribution function Fi(t) based on the position information at time t.

1We are interested in second order statistics because of the deﬁnition of emittance. A detailed discus-

sion will be presented in the Sec. 3.2

35

Therefore, the dynamics of these three statistics can be expressed as:

dσ2
i
dt
dσ2
pi
dt
dσi,pi
dt

2
m

(cid:104)ipi(cid:105) =

(cid:39) 2
σi,pi
m
(cid:39) 2(cid:104)piFi(cid:105) = 2σpi,Fi
(cid:39) (cid:104)iFi +

1
m

i(cid:105) = σi,Fi
p2

(3.3a)

(3.3b)

(3.3c)

+

1
m

σ2
pi

We also have another set of three statistics speciﬁcally for SSA, which is equivalent to

Eq. 3.1, as follows:

i = σ2
σ2
i

γi = σi,pi

ηi = σ2
pi

− σ2
i,pi
σ2
i

(3.4a)

(3.4b)

(3.4c)

where ηi is the variance of local momentum ﬂuctuations, the variance of the momentum

after subtracting out the portion that is linearly correlated with position as shown below.

A detailed deﬁnition2 of these three statistics is illustrated in Fig. 3.1. In the longitudinal

direction, σ2

z is the bunch spatial variance; γz is the correlation between the longitudi-

nal position and momentum, and ηz is the variance of local momentum ﬂuctuations in

the longitudinal direction. Speciﬁcally, for each individual electron [j], its longitudinal
momentum pz[j] can be expressed as a linear function (with slope Cz) of its longitudinal

2More underlying signiﬁcance of these three statistics can be found in Sec. 3.3.1.

36

Figure 3.1: Schematic representation for the three bunch parameters in SSA

position z[j] plus a random variable δz[j] as follows:

pz[j] = Cz · z[j] + δz[j] =

γz
σ2
z

z[j] + δz[j]

(3.5)

Here δz serves as the error term, with variance ηz = σ2

δz. In addition, δz is assumed to

be uncorrelated with the longitudinal position, i.e. σδz,z = 0. With this set of notation,

SSA essentially divides the variance of momentum into two parts: one part comes from

the momentum which is linearly dependent on position of each particle and the other part

is associated with the momentum that deviates from the linear part. The corresponding

37

dynamics for SSA statistics can be derived from Eq. 3.3 as follows:

γi

(cid:32)

dσ2
i
dt
dγi
dt

=

=

2
m
1
m

ηi +

dηi
dt

= −2γiηi
mσ2
i

+

γ2
i
σ2
i
2
σ2
i

(cid:33)
(cid:16)

+ σi,Fi

σ2
i σpi,Fi

− σi,pi

σi,Fi

(cid:17)

(3.6a)

(3.6b)

(3.6c)

The necessary initial conditions can be obtained by gathering the second-order statistics

of the bunch density distribution at the beginning of the evolution, namely the three

bunch statistics σ2

i , γi, and ηi can be derived from the three direct second-order statistics

of σ2

i , σi,pi

, and σ2
pi

.

Coupled with proper force term, the SSA prediction of bunch evolution can be obtained

by numerically propagating the SSA dynamics from the initial conditions. Therefore,

the key of SSA — that the bunch should evolve self-similarly — comes from the force

distribution Fi. In other words, SSA assumes that the uniform ellipsoidal bunch will stay

uniform while the density proﬁle of a Gaussian bunch will stay as a Gaussian function.

The self-similar evolution assumption simpliﬁes the analysis by setting the force expression

term Fi as the same form of function throughout the evolution with one (or more) time-

dependent parameters in the charge density distribution.

For a uniform spheroidal bunch, as shown in Eq. 2.17, the electrostatic ﬁeld is linearly

dependent on the spatial position. In other words, the force term can be written as:

Fi =

ρce
2ε0

· ξi(α) · i = KF i(α) · i

(3.7)

38

where i = T, z for the transverse3 and longitudinal directions and α is the transient aspect

ratio of the spheroid at time t. The corresponding second-order statistics related to Fi

can then be simpliﬁed as:

σi,Fi

σpi,Fi

= KF i(α) · σ2
= KF i(α) · σi,pi

i

which also leads to the vanishing of the second term in Eq. 3.6c as:

σ2
i σpi,Fi

− σi,pi

σi,Fi

= 0

(3.8a)

(3.8b)

(3.9)

Therefore, the resulting SSA prediction of the self-similar evolution of a uniform spheroidal

bunch yields:

dσ2
i
dt
dγi
dt

dηi
dt

=

2
m
1
m

γi

(cid:32)

=

γ2
i
σ2
i
= −2γiηi
mσ2
i

(cid:33)

+ ηi

+ KF i(α) · σ2

i

(3.10a)

(3.10b)

(3.10c)

One advantage of SSA is that the two-step process for dealing with the singularity at

longitudinal crossover in MFT is not necessary here in the SSA model as the statistics

themselves are second-order moments so that there is no need to ﬂip the sign of spatial

position and momentum in the cross-over regime.

3As we take advantage of the rotational symmetry of spheroidal ellipsoids, the subscript T represents

the two transverse directions, which includes both ˆx and ˆy.

39

Although the above derivations have been performed for Coulomb interactions, we

would like to stress that the same conclusions can be drawn for any interaction that leads

to linear dependence between force and position, such as gravitation [31]. Additionally, the

generalization to any general ellipsoid is achieved by simply using three degrees of freedom

with i = X, Y, Z and corresponding geometry coeﬃcients (ξx, ξy, ξz) with functions of the

ratio between three axes (σx : σy : σz).

3.2 Bunch evolution with conserved emittance

One of the major advantages of the SSA model (Eq. 3.10) is its inclusion of a conserved

emittance, while MFT always assumes a zero-emittance bunch. In this analysis, we deﬁne

the “root-mean-square” (rms) emittance εi in each direction as:

ε2
i = σ2

i σ2
pi

− σ2
i,pi

= σ2

i ηi

(3.11)

This suggests that the value of emittance in each direction is proportional to the area of

the contouring ellipse of the bunch in the 2D phase space regarding that direction. In

other words, the rms emittance is not equal to the exact value of the occupied phase-

space volume, but rather a statistical estimate of bunch phase-space volume based on the

second-order statistics of the bunch. Although some literature conﬂates rms emittance

with phase-space volume, these two concepts, while related, are not equal to each other.

According to [43], rms emittance can grow while the actual phase-space volume stays

constant as explained by Liouville’s theorem.

To facilitate the analysis, we base our derivation using the square of emittance ε2
i

40

which is the product of σ2

i and ηi. Thus, the evolution of squared emittance is:

d
dt

ε2
i = ηi

d
dt

i + σ2
σ2
i

d
dt

ηi

(3.12)

Although, we use the squared emittance ε2

i for analytical derivations and discussions, we

present results in emittance εi to ease the comparisons between model prediction and

actual experiments.

As we look at the emittance evolution for uniform bunches, it is conserved within

each degree of freedom according to Eq. 3.10, which seems to be the immediate result

of the linear ﬁeld/force assumption of Eq. 3.7. However, as we elaborate later in the

next section, it is the linear chirp (self-similar evolution) assumption that ensures the

conservation of emittance of the bunch.

Additionally, if we take η = 0, (i.e. zero emittance), then the SSA model will reproduce

the MFT, as shown in Fig. 3.2. The linearity/self-similar evolution assumption is the

reason why SSA, as a statistical method, can reproduce the result of MFT. The self-

similar evolution assumption ensures the linearity across the bunch for SSA, as the linear

momentum-position relationship in Eq. 3.5.

The emittance eﬀect on the SSA dynamics can be explained by studying the evolution

of the linear chirp Ci as we have:

Ci =

d
dt

d
dt

(cid:33)

(cid:32)

γi
σ2
i

= KF i(α) +

m ·(cid:0)σ2

ε2
i

i

(cid:1)2

(3.13)

During the focusing process, the emittance term (last term in Eq. 3.13), will drive the

negative chirp to approach zero faster. Therefore, the waist of the focusing process ap-

41

c , (b)1.0ω∗

Figure 3.2: Longitudinal width evolution of prolate spheroids with (α0 = 10/75) focused
by diﬀerent initial chirps: (a) 0.7ω∗
c . In each ﬁgure, red solid line
represents the prediction from MFT, dotted lines represent SSA with diﬀerent emittance
ranging from 0 to 0.01µm or mm · mrads. The SSA with zero emittance and the MFT
are in perfect agreement. Notice that (1) emittance move the waist larger and earlier
comparing SSA and MFT, (2) comparing diﬀerent chirps in SSA, the bunch evolution
driven by critical chirp(b) shows higher sensitivity for emittance.

c , (c) 1.5ω∗

42

0.50.60.70.80.91.0(a)0.00.20.40.60.81.0z(t)/z(0)(b)0100200300400500600700t (ps)0.00.20.40.60.81.0(c)MFTSSA(z=10.00nm)SSA(z=3.00nm)SSA(z=1.00nm)SSA(z=0.30nm)SSA(z=0.10nm)SSA(z=0.00nm)pears earlier than that predicted by MFT, resulting in a larger minimum width as well.

Moreover, bigger emittance will lead to larger diﬀerences between predictions from SSA

and MFT. These two emittance eﬀects are also shown in Fig. 3.2.

3.3 Linear chirp assumption conserves emittance

To further illustrate the fact that the linear chirp assumption conserves emittance in the

SSA model, we provide the following derivation.

Although the electron bunches can have various density proﬁles, such as uniform or

parabolic distributions, we start with an electron bunch with Gaussian density proﬁle as

an example. We rely on this example because the expression of the Gaussian density

proﬁle can have three explicit bunch statistics, the same set as in SSA, which makes

it more advantageous for our analysis. The expression for the Gaussian density proﬁle

follows this general form:

f (r, p; t) = C(t) exp [Σi (Gi(i) + Hi(pi, i))]

(3.14)

where i = x, y, z. Gi is the function for the spatial distribution in each direction, which

can be written as:

Gi(i) =

a0 + a2i2 + a4i4 + ... + a2n−2i2n−2 − a2ni2n

2σ2
i

(3.15)

43

Meanwhile, Hi is the function for the momentum distribution, which can be written as:

(cid:2)pi −(cid:0)γi/σ2

(cid:1) i(cid:3)2

Hi(pi, i) = −

i

2ηi

(3.16)

where (σ2

i , γi, ηi) are the time-dependent statistical descriptions of the bunch in each

direction we previously deﬁned in SSA. The linear chirp of the bunch can still be described

as Eq. 3.5, which leads to the self-similar evolution. The corresponding time derivative

of squared emittance reads:

d
dt

ε2
z = ηz

dσ2
z
dt

+ σ2
z

dηz
dt

= 2

(cid:16)

z σpz,Fz − σz,pz σz,Fz
σ2

(cid:17)

(3.17)

3.3.1 The Linear chirp assumption in detail

In the proceeding section, I concluded that uniform ellipsoidal bunches conserve emittance

after inserting the corresponding linear force from its self-ﬁeld into the SSA model. Here

I show that the linear chirp assumption also leads to conserved emittance, which is not

restricted to the uniform bunch density proﬁle as a linear force assumption. The linear

chirp assumption focuses on the statistical relationship between the spatial position and

momentum for all the electrons of the bunch. Speciﬁcally, the linear chirp assumption

assumes that the momentum of each individual electron j can be divided into two parts:

the local average momentum based on the spatial position of electron j and the ﬂuctuation

part for electron j, just as we described in Eq. 3.5:

pz[j] = Cz · z[j] + δz[j] =

γz
σ2
z

z[j] + δz[j]

44

In the linear chirp assumption, the local average momentum is proportional to its spatial

coordinates. That is, for longitudinal direction, the slope of this linear dependence (be-
tween local average momentum and spatial coordinate) Cz is decided by the two bunch
statistics as γz
σ2
z

. Then, δz[j] is denoted as the diﬀerence between electron j’s momentum

and the projected local average momentum that is linearly based on its spatial coordinate,

as shown in Eq. 3.5. We refer to δz as the local momentum ﬂuctuation in the longitudinal

direction with respect to the corresponding linear local average momentum. The variance

of this local momentum ﬂuctuation is then equal to ηz as we have:

=

i+δz

γ2
z
σ4
z

σ2
z + 2

γz
σ2
z

σz,δz + σ2
δz

(3.18)

+ ηz

σ2
pz = σ2
γz
σ2
z
γ2
z
σ2
z

=

with the assumption that δz and z are uncorrelated, i.e. σz,δz = 0.

This assumption about δz and z being uncorrelated is the center piece for the linear

chirp assumption. Because the linear local average momentum can always be projected

to the bunch density proﬁle from the corresponding bunch statistics, the linear chirp

assumption is essentially assuming the local momentum ﬂuctuation with respect to that

linear average momentum is uncorrelated with the spatial position. In other words, the

momentum of electrons in bunches that do not satisfy that linear chirp assumption can

still be expressed in Eq. 3.5, but the corresponding local momentum ﬂuctuation (δi) will

then be correlated with the spatial position. Therefore, the linear chirp assumption is

actually focusing on the local momentum ﬂuctuation (with respect to that linear chirp)

being uncorrelated with the spatial position.

45

One possible expression of H in a Gaussian density proﬁle that satisﬁes the linear

chirp assumption as Eq. 3.16:

Hz(pz, z) = − 1
2ηz

(cid:20)

pz − γz
σ2
z

(cid:21)2

z

which assumes that for any point inside the ellipsoidal bunch, the density distribution

in momentum space with respect to the corresponding local average momentum has a

Gaussian proﬁle. The two bunch statistics involving pz become:

σpz,Fz = σ γz
σ2
z

z+δz,Fz

=

γz
σ2
z

σz,Fz + σδz,Fz

σz,pz = σz,

γz
σz z+δz

=

γz
σz

σ2
z + σz,δz

(3.19a)

(3.19b)

(3.20)

(cid:19)

(cid:21)

σz,Fz

(cid:18) γz

σz

−

(cid:19)
(cid:17)

σ2
z + σz,δz

Accordingly, we can rewrite the evolution of emittance as:

(cid:18) γz

(cid:20)
(cid:16)

dε2
z
dt

= 2

= 2

σ2
z

σz

σz,Fz + σδz,Fz
z σδz,Fz − σz,δz σz,Fz
σ2

The second term equals to zero as σz,δz = 0 from the linear chirp assumption. The ﬁrst

term also vanishes if the bunch density proﬁle is symmetric with respect to the spatial

position and ﬂuctuation of momentum. Speciﬁcally, we can have:

˜f (z, δz) = ˜f (z,−δz) = ˜f (−z, δz)

(3.21)

46

with ˜f being the shifted density distribution of the original one as:

˜f (z, δz) = f (z,

γz
σ2
z

z + δz) = f (z, pz(z, δz))

(3.22)

for a symmetric bunch density proﬁle. The immediate consequence of a symmetric density

distribution is that the magnitude of force due to the self-ﬁeld can be expressed as an odd

function of the spatial coordinates, i.e.

Fz(x, y, z) = −Fz(x, y,−z)

(3.23)

Then, the following covariance term related to the force distribution vanishes as:

(cid:90) (cid:90)

σδz,Fz =

δzFz(z) ˜f (z, δz; t)drdp

(3.24a)

Speciﬁcally, this covariance term vanishes when we pair terms like δzFz(z) ˜f (z, δz; t) and
δzFz(−z) ˜f (−z, δz; t) together, as these cancel each other out. As both terms in Eq. 3.20

equal to zero for a symmetric density proﬁle, the emittance is conserved throughout the

evolution, as:

dε2
z
dt

= 0

(3.25)

In addition, any symmetric density proﬁle that satisﬁes Eq. 3.21 can lead to conserved

emittance, not limited to the general Gaussian density proﬁle mentioned at the beginning

of this section.

As demonstrated throughout these derivations, it becomes clear that the exact form

of the electron bunches’ density proﬁles do not aﬀect emittance conservation [38]. As

47

long as Eq. 3.21 is satisﬁed, the emittance is conserved. In other words, if the ﬂuctuation

of momentum δi of each individual electron is uncorrelated with its spatial position, the

ellipsoidal symmetric bunch conserves emittance throughout its evolution.

3.3.2 A Comparison Between the SSA model and the Analytic

Gaussian model

To further validate the argument that emittance conservation does not depend on the

exact form of the spatial portion of the bunch’s density proﬁle, we overview the Analytic

Gaussian (AG) model — a semi-statistical mean-ﬁeld approach for bunch evolution with

Gaussian density proﬁles.

Michalik and Sipe [34] introduced the AG model to study the evolution of bunch

statistics with Gaussian proﬁles in both spatial and momentum space. Their key assump-

tion is that bunches retain their Gaussian density proﬁles with the three time-dependent

parameters in each direction. We generally base the following derivation on the original

AG paper, however, we substitute in the SSA notation for consistency throughout this

analysis.

We start with the Gaussian density proﬁle used in the AG model:

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

(3.26)

48

with

Γ(r, p; t) =

+

x2 + y2

2σ2
T
z2
2σ2
z

+

+
[pz − (γz/σ2

z )z]2

2ηz

[px − (γT /σ2

T )x]2 + [py − (γT /σ2

T )y]2

2ηT

(3.27)

This expression has three main assumptions. First, this expression assumes rotational

symmetry in ˆx and ˆy direction. The subscript T then represents the two transverse

directions, which includes both ˆx and ˆy. The longitudinal direction is denoted with

subscript ˆz. Second, the three bunch parameters in the AG model,4 (σ2

z , γz, ηz) are the

same as those used in the SSA model. This expression means the bunch has linear chirps

in each direction and the momentum follows a relationship similar to Eq. 3.5, which is

illustrated in Fig. 3.1 for the longitudinal direction. In other words, the Gaussian bunch

in the AG model has linear chirp in each direction, with the corresponding momentum

spread is uncorrelated with the spatial position. Therefore, the AG model is essentially

modeling the longitudinal focusing of a spheroidal Gaussian bunch by calculating the self-

similar evolution of a Gaussian bunch with the same initial bunch statistics. Third, the

time-dependent normalization factor C(t) comes from the conservation of total charge,

which can be expressed as:

(cid:90)

f (r, p; t)drdp = N ⇒ C(t) =

(cid:32)

N

(2π)3

1
σ2
T η2
T σzηz

(cid:33)1/2

(3.28)

That is, this density proﬁle is a function of the three parameters in each direction (T and

z). Therefore, the evolution of Gaussian bunches in the AG model can be generated from

4We are illustrating the parameters for longitudinal direction. Similar notations apply to transverse

direction as well.

49

the trajectories of the three parameters of the bunch in each direction.

By assuming the Gaussian bunch follows a self-similar evolution, the AG model allows

researchers to model bunch evolution with Coulomb interactions in two main ways. First,

the impact of the Coulomb interactions can be modeled by evaluating the corresponding

Gaussian integrals. Speciﬁcally, the bunch evolution in each direction can be expressed

by the trajectories of the three parameters in Eq. 3.27. The resulting force terms in

the time-derivatives of the bunch parameters are expressed as functions of the transient

aspect ratio (α(t) = σT /σz). Second, because the bunch retains a Gaussian proﬁle, the

AG model captures the Coulomb eﬀect throughout the evolution of the electron bunch.

Speciﬁcally, this can be expressed using the same force term with the time-dependent

aspect ratio. These two advantages of the AG model make the theoretical expression of

the bunch evolution with Coulomb interaction feasible assuming a self-similar Gaussian

bunch.

The pathway to the force term is a two-step process. First, we calculate the three

parameters of the Gaussian density proﬁle f (r, p; t) to calculate the Coulomb eﬀect on

the bunch emittance evolution. Second, we integrate over the “ﬂow” term and “force”

term from the time derivative of f to eventually solve for the time derivative of these

three parameters.

50

Speciﬁcally, for the longitudinal direction,5 we have:

σ2
z =

γz =

γ2
z
σ2
z

=

1
N
1
N
1
N

ηz +

z2f (r, p; t)drdp = X55
zpzf (r, p; t)drdp = X56
zf (r, p; t)drdp = X66
p2

(3.29a)

(3.29b)

(3.29c)

(cid:90)
(cid:90)
(cid:90)

(cid:90)

Here we use Xij to represent the corresponding second moment of the Gaussian proﬁle

f (r, p; t) for simplicity. That is:

Xij =

1
N

(cid:90)

uiujf (r, p; t)drdp

(3.30)

with (u1, u2, u3, u4, u5, u6) = (x, px, y, py, z, pz). Conveniently, we can have the time
derivative of Xij as:

Xij =

∂
∂t

1
N

uiuj

∂f (r, p; t)

∂t

drdp

Therefore, the time derivatives of the three parameters can be written as:

(cid:90)
(cid:90)

1
N
1
N

∂
∂t
∂
∂t
∂
∂t

σ2
z =

z2 ∂f (r, p; t)

∂t

drdp =

γz =

zpz

ηz = −2γz
σ2
z

(cid:18) γz

∂
∂t

(cid:19)2 ∂

=

σ2
z

∂t

∂f (r, p; t)

γz +

∂t

γ2

(cid:1)2
z(cid:0)σ2
X55 − 2

z

X55
(cid:90)
X56

∂
∂t
∂
∂t
1
N
X56 +

drdp =

∂
∂t

σ2
z +

γz
σ2
z

∂
∂t

X66

∂
∂t

p2
z

∂f (r, p; t)

∂t

drdp

5Similar results can be applied to the transverse direction for both ˆx and ˆy.

51

(3.31)

(3.32a)

(3.32b)

(3.32c)

According to the Vlasov equation, we can rewrite ∂f /∂t as:

∂f (r, p; t)

∂t

= −∂f (r, p; t)

(cid:16)− p

∂r

∂r
∂t

(cid:17) ∂f (r, p; t)

− ∂f (r, p; t)

∂p
∂t
− ∂f (r, p; t)

∂p

=

m

∂r

∂p

F (r)

(3.33)

The integral of the ﬁrst term is called the “ﬂow” term, Kf low

ij

, as it comes from funda-

mental kinematics. The integral of the second term is called the “force” term, Kf orce

ij

.

The force term represents the eﬀect of the Coulomb interaction on the bunch evolution.

(cid:90)
(cid:90)

Accordingly, we rewrite Eq. 3.31 as:

Xij =

∂
∂t

=

∂t

uiuj

(cid:16)− p

1
N
1
N
m
ij + Kf orce

uiuj

ij

= Kf low

∂f (r, p; t)

drdp

(cid:17) ∂f (r, p; t)

∂r

(cid:90)

uiuj

drdp +

1
N

(cid:18)

−∂f (r, p; t)

∂p

(cid:19)

F (r)drdp

(3.34)

We start with the ﬂow term. The spatial derivative of the Gaussian density proﬁle f

is:

∂f (r, p; t)

∂z

=

=

(cid:20) 1
(cid:20)(cid:18) 1

σ2
z

(cid:18)
(cid:19)

(cid:19)(cid:21)
(cid:21)

z − γz
σ2
z ηz
γ2
z
σ4
z ηz

+

σ2
z

z

pz − γz
σ2
z
z − γz
σ2
z ηz

f (r, p; t)

pz

f (r, p; t)

(3.35)

52

Therefore, the corresponding ﬂow terms are:

(cid:90) (cid:90)
(cid:90) (cid:90)
(cid:90) (cid:90)

z2pz

zp2
z

σ2
z

(cid:20)(cid:18) 1
(cid:20)(cid:18) 1
(cid:20)(cid:18) 1

σ2
z

Kf low

55 =

Kf low

56 =

1

N m

1

N m

Kf low

66 =

1

N m

+

γ2
z
σ4
z ηz
γ2
z
σ4
z ηz

+

z − γz
σ2
z ηz
z − γz
σ2
z ηz

(cid:19)
(cid:19)
(cid:19)

pz

(cid:21)
(cid:21)
(cid:21)

f (r, p; t)drdp =

(3.36a)

pz

f (r, p; t)drdp =

ηz +

γz

2
m

(cid:18)

1
m

(cid:19)

γ2
z
σ2
z

(3.36b)

p3
z

γ2
z
σ4
z ηz

+

σ2
z

z − γz
σ2
z ηz

pz

f (r, p; t)drdp = 0

(3.36c)

If we assume F (r) = 0, then the time derivatives of the three parameters (Eq. 3.32), with

the eﬀect of the ﬂow term alone, models the evolution of a non-interacting bunch with a

Gaussian density proﬁle as Eq. 3.27. That is, the non-interacting bunch evolution could

be described by inserting Eq. 3.36 in the Eq. 3.32 as:

dσ2
z
dt
dγz
dt
dηz
dt

γz

(cid:18)

ηz +

=

=

2
m
1
m

(cid:19)

γ2
z
σ2
z

= −2γzηz
mσ2
z

(3.37a)

(3.37b)

(3.37c)

Since the three bunch statistics are only functions of time, we replace the partial deriva-

tives with the total derivatives here to obtain the ordinary diﬀerential equation for the

bunch statistics evolution. Unsurprisingly, the non-interacting bunch conserves emittance

throughout the evolution when we insert Eq. 3.37 into Eq. 3.12.

For the purpose of measuring the Coulomb interaction eﬀect on bunch evolution,

especially on emittance conservation, we insert the corresponding force expression F (r)

derived from the expression of Coulomb potential, which we derived earlier (Eq. 2.47).

53

Consequently, the time derivative of the three parameters (Eq. 3.32) yields the collective

eﬀect of a Gaussian bunch through Kf orce

ij

.

Speciﬁcally, the ﬁeld distribution can be derived from Eq. 2.47 as:

(cid:90) ∞

Q

4π3/2ε0

0

2z

c2 + l

·

Ez(x, y, z) =

Employing the following substitutions:

(cid:20)
(cid:112)(a2 + l)(b2 + l)(c2 + l)

− z2
c2+l

− y2
b2+l

− x2
a2+l

exp

(cid:21)

dl

(3.38)

(3.39a)

(3.39b)

(3.39c)

(3.39d)

Qe

2π3/2ε0

K =

R(l) =

λ(x, y, z; l) =

Ic(x, y, z) =

(cid:113)
(cid:90) ∞

x2

(a2 + l)(b2 + l)(c2 + l)

y2

z2

a2 + l

+

b2 + l

+
c2 + l
exp [−λ(x, y, z; l)]

0

(c2 + l)R(l)

dl

we can derive a concise expression for the Coulomb force as:

Fz(x, y, z) = z · K · Ic(x, y, z)

(3.40)

Then, knowing that the derivative of f with respect to pz is

∂
∂pz

f (r, p; t) =

1
ηz

pz − γz
σz

z

f (r, p; t)

(3.41)

(cid:19)

(cid:18)

54

the force term in the longitudinal direction can be expressed as:

(cid:90) (cid:90)
(cid:90) (cid:90)
(cid:90) (cid:90)

(cid:18)

(cid:19)

z3

z2pz

zp2
z

(cid:19)
(cid:19)

z

z

(cid:18)
(cid:18)

pz − γz
σ2
z
pz − γz
σ2
z
pz − γz
σ2
z

z

Kf orce

55

Kf orce

56

Kf orce

66

= − 1
N ηz
= − 1
N ηz
= − 1
N ηz

· Ic(x, y, z) · f (r, p; t)drdp

· Ic(x, y, z) · f (r, p; t)drdp

· Ic(x, y, z) · f (r, p; t)drdp

(3.42a)

(3.42b)

(3.42c)

(3.43)

For the emittance evolution, we can insert Eq. 3.32 into Eq. 3.12 as:

dε2
z
dt

= ηz

(cid:18)

dσ2
z
dt

=

ηz +

+ σ2
z

(cid:19)

γ2
z
σ2
z

dηz
dt
Kf orce

55

+ σ2
z

(cid:18)

66 − 2γz
Kf orce
σ2
z

(cid:19)

Kf orce

56

Thus the emittance evolution with Eq. 3.42:

(cid:90) (cid:90) (cid:34)

(cid:18)

(cid:19)

(cid:18)

dε2
z
dt

= − 1
N

z3

pz − γz
σ2
z

z

+

σ2
z
ηz

z

pz − γz
σ2
z

(cid:19)3(cid:35)

z

· Ic(x, y, z) · f (r, p; t)drdp

We can simplify this expression by a substitution of variable as:

δz = pz − γz
σ2
z

z

with the corresponding Jacobian:

J = det

 1

− γz
σ2
z

 = 1

0

1

(3.44)

(3.45)

(3.46)

Accordingly, the phase-space area of z − pz are the same as that of the corresponding

55

z − δz phase-space. The density distribution function f (r, p; t) shifts to ˜f as shown in

Eq. 3.22. We can rewrite the time derivative of ε2

z as:

(cid:90) (cid:90) (cid:20)

dε2
z
dt

= − 1
N

(cid:21)

z3δz +

σ2
z
ηz

zδ3
z

· Ic(x, y, z) · ˜f (r, δr; t)drdδr

(3.47)

Since Ic(x, y, z) and ˜f are both even functions, this integral equals zero as we properly

pair the elements as expressed below:

(z, δz)

and

(−z, δz)

(3.48)

Therefore, emittance is conserved during the evolution of a Gaussian bunch according to

the AG model, even with a non-linear self-ﬁeld/force. This is consistent with what we

argued in the proceeding section, that the linear chirp is the key assumption for conserved

emittance, not the linear force.

We continue the derivation for the Coulomb eﬀect on bunch evolution in the AG
model. To further simplify the expression, we introduce the “bar” notation (Mz, Mz as a

place-holder for the corresponding quantities, which will be evaluated for the longitudinal

direction) as:

Mz = − 1
N ηz

(cid:90) (cid:90)

Mz · Ic(x, y, z) · ˜f (r, δr; t)drdδr

(3.49)

56

Thus, we can rewrite the force term in Eq. 3.42 as

Kf orce

55

Kf orce

56

Kf orce

66

(cid:18)
(cid:18)

= z3δz = 0

= z2δz

δz +

γz
σ2
z

= zδz

δz +

γz
σ2
z

z

(cid:19)
(cid:19)2

z

= z2δ2

z +

= zδ3

z + 2

γz
σ2
z
γz
σ2
z

z3δz = z2δ2
z

z2δ2

z +

(cid:18) γz

σ2
z

(cid:19)2

(3.50a)

(3.50b)

z3δz = 2

γz
σ2
z

z2δ2

z =

2γz
σ2
z

Kf orce

56

(3.50c)

where the integrals of z3δz and zδ3

z vanish due to the symmetry of Ic(x, y, z) and ˜f (r, δr; t)

(i.e. f (r, p; t)). As we can have a closed form z2δ2

z and similarly for x2δ2

x and y2δ2
y:

Kf orce

56

=

Kf orce

12

= Kf orce

34

=

1

4π0

1

4π0

N e2
√
6
π
N e2
√
6
π

1
σz
1
σT

Lz

LT

(cid:19)
(cid:18) 1
(cid:18) 1
(cid:19)

α

α

(3.51a)

(3.51b)

(3.52a)

(3.52b)

where Li(p) is:

Lz(p) =

LT (p) =

and L(p) as:

(cid:90) π

dθ

1 + p sin θ

0

L(p) =

1
2



=

(cid:20)

3p2
p2 − 1
3
2

L(p) +

[pL(p) − 1]

p2L(p) − p

1 − p2

(cid:21)

arcsin

1−p2

0 ≤ p ≤ 1 (oblate)

,

ln(p+

p2−1)

1 ≤ p (prolate)

,

(3.53)

(cid:113)
(cid:113)

(cid:113)
(cid:113)

1−p2

p2−1

57

Once the ﬂow term and force term are known, we can write the dynamics expression

for the self-similar evolution of a Gaussian bunch with the three parameters:

(cid:18) 1

(cid:19)

α

+

1

4π0

N e2
√
6
π

1
σi

Li

dσ2
i
dt
dγi
dt

=

=

2
m
1
m

γi

(cid:32)

ηi +

(cid:33)

γ2
i
σ2
i

dηi
dt

= −2γiηi
mσ2
i

(3.54a)

(3.54b)

(3.54c)

3.3.3 Discussion about AG

The essence of the AG model is the evolution of a self-similar Gaussian bunch, where

the chirp is linear and the local momentum ﬂuctuation is uncorrelated with the electron’s

spatial position.

It uses the initial bunch statistics as the initial values of the three

time-dependent parameters of the Gaussian proﬁle. Therefore, the credibility of AG’s

prediction comes down to the similarity between the actual charge density proﬁle and the

Gaussian proﬁle that shares the same bunch statistics. This is the reason behind the good

performance for the AG model, despite the fact that a Gaussian charge distribution does

not evolve self-similarly [17, 36] due to the non-linear force.

To this extent, the AG model can be considered as an intermediate approach, which

shifts the focus from the exact density proﬁle towards bunch statistics. In other words,

the advantage of the AG model is that it starts to focus more on the evolution of bunch

statistics. However, the AG model is limited in that is assumes that bunches have and

retain a Gaussian form.

58

3.4 Discussion

We presented the SSA model, a statistical method describing the second order moments to

model the evolution of an ellipsoidal electron bunch and explained how it can reproduce

the MFT by taking a zero momentum spread (i.e. zero emittance). We also discussed

the linear chirp assumption, explaining how it is the key assumption to yield emittance

conservation during the bunch evolution. In addition, we also shed light on the statistical

nature of emittance, noting how the force impacts from the uniform density proﬁle and

Gaussian proﬁle are similar to each other. In the following subsections, we expand our

discussion upon all of these points.

3.4.1 linear chirp and self-similar evolution assumptions

The derivation outlined in the preceding section of this chapter showed that the key

assumption for conserving emittance is the linear chirp assumption. In particular, the

linear local average momentum can always be projected on the bunch distribution function

by inserting the corresponding statistics into Eq. 3.5. The ﬂuctuation of momentum (δi)

is then the crucial quantity that diﬀerentiates bunches from conserving emittance or not.

Just as we mentioned in the discussion about the MFT, if δi is not symmetric, or not

uncorrelated with the spatial position, we can still use the SSA model or the AG model

to estimate the evolution. However, these model estimations will deviate from the true

evolution. In the next chapter, we explore the magnitude of this discrepancy.

We refer to the SSA model as the self-similar analytic model to diﬀerentiate it from

the AG model. Speciﬁcally, through this naming, we emphasize the self-similar evolution

assumption rather than the Gaussian density proﬁle.

In this analysis, the self-similar

59

bunch evolution is also considered satisfying the linear chirp assumption throughout the

evolution as this ensures that the bunch with a linear chirp will evolve self-similarly.6 In

reality, the Gaussian bunch generally does not evolve self-similarly, as pointed out by the

literature [17,36] as is also demonstrated in the simulation in the following Sec. 4. Indeed,

only uniform bunches with zero emittance evolve truly self-similarly.

3.4.2 Kinetic energy transfer between degrees of freedom

The kinetic energy statistics of the bunch in each degree of freedom is related to the three

statistics through:

KEi =

N
2m

σ2
pi

=

N
2m

ηi +

γ2
i
σ2
i

(cid:32)

(cid:33)

(3.55)

The corresponding result of a crossover case is shown in Fig. 3.3 The kinetic energy

evolution of a focusing ellipsoidal electron bunch can be written as:

d
dt

KEi =

γiKF i(α)

N
m

(3.56)

This evolution is consistent with the fundamental idea that focusing an electron bunch

against Coulomb repulsion transfers kinetic energy in the focusing degree of freedom

into potential energy while expansion works in the other directions. For example, when

the electron bunch is focused in the longitudinal direction, the kinetic energy associated

with the longitudinal direction will transfer to the potential energy. Therefore, if one

degree of freedom is focusing while the other one is expanding, it would appear to be

6There will be edge eﬀect due to the ﬂuctuation of momentum that create density tails on the bunch
edges, which means the bunch evolution is not strictly self-similar. We present a detailed discussion in
the following section.

60

Figure 3.3: Kinetic energy evolution for cross-over case ( −1.5ω∗
c in the longitudinal
direction, corresponding to panel (c) in Fig. 3.2), with solid lines for MFT and dotted
lines for SSA with diﬀerent value of emittance (circle for KEz and triangle for KEr).
The sudden change of direction for MFT in longitudinal kinetic energy comes from the
sign ﬂip of the chirp as discussed in Sec. 2.

kinetic energy transfers between two directions. This energy transfer is mediated by the

potential energy as the kinetic energy in the focusing direction goes into the potential

energy and the potential energy converts to the kinetic energy in the expanding direction.

Therefore, we can control this energy transfer process by tuning the focusing procedure.

3.4.3 Statistical discussion and the K-V envelope equations

In this subsection, we want to discuss the statistical essense of the SSA model, namely that

the SSA model describes the evolution of second-order moments of the bunch distribution.

We compare the force term for both a uniform bunch (Eq. 3.10) and a Gaussian bunch

61

050100150200250t (ps)050100150200KE(Hartree)KEz[MFT]KEr[MFT]KE[SSA](z=10.00nm)KE[SSA](z=3.00nm)KE[SSA](z=1.00nm)KE[SSA](z=0.00nm)(Eq. 3.54), as the force term for a uniform spheroidal electron bunch can also be evaluated

analytically.

Let us take the longitudinal direction of a uniform oblate (a = b > c) spheroidal bunch

as an example. For a uniform spheroidal bunch, we can have the following translations

between parameters and bunch statistics:

√
√

5σT

5σz

3N
4πa2c

a =

c =

n =

(3.57a)

(3.57b)

(3.57c)

The linear coeﬃcient KF z(α) in the force terms in the longitudinal direction from Eq.
3.10 can be rewritten as:7

a2c

0
N e2
3
c3
2
N e2
√
5σ3
5
z

(cid:90) ∞
(cid:90) ∞
(cid:90) ∞
 2c2
(cid:18) 1
(cid:19)

N e2
√
5σ3
5
z

0

0
3
2

3
2

N e2
√
5σ3
5
z

Lz

α

KF z(α) =

=

=

=

=

ne2
2ε0
1

4πε0

1

4πε0

1

4πε0

1

4πε0

dl

(c2 + l)3/2(a2 + l)

(1 + l/c2)3/2(a2/c2 + l/c2)

d(l/c2)

dν



(1 + ν)3/2(a2/c2 + ν)

c2 − a2

c√
a2 − c2

arcsin

(cid:115)



1 − c2

a2 − 1

(3.58)

Similar results apply to the transverse direction and both directions for prolate spheroids

as well. We can then rewrite the SSA evolution equations for a uniform spheroidal bunch

7Credits to Brandon Zerbe for the evaluation of the corresponding integrals.

62

as:

(cid:18) 1

(cid:19)

α

+

1

4π0

N e2
√
5
5

1
σi

Li

γi

(cid:32)

dσ2
i
dt
dγi
dt

=

=

2
m
1
m

ηi +

γ2
i
σ2
i

(cid:33)

dηi
dt

= −2γiηi
mσ2
i

(3.59a)

(3.59b)

(3.59c)

with i = T, z for transverse and longitudinal direction, respectively. The impact of the

collective Coulomb interaction from a uniform density proﬁle on bunch evolution is strik-

ingly similar to that from a Gaussian bunch (Eq. 3.54). The diﬀerence between the two

density proﬁles can be expressed by the ratio of the two force terms as:

KGaussian
KU nif orm

=

√
√
5
6

5
π

≈ 1.0513

(3.60)

Therefore, the diﬀerence between the two density proﬁles is only a little over 5%. In other

words, according to the SSA model, the evolution of bunch statistics for a uniform bunch

with N electrons is identical to that of a Gaussian bunch with the same initial bunch

statistics and total number of electrons equals to 0.95N .

This similarity comes from two key factors: (1) σi,Fi

describes the linear part of the

self-ﬁeld with respect to the corresponding spatial position, which is the only presence

of the Coulomb interaction for bunch evolution in the SSA model with the linear chirp

assumption (2) the mean-ﬁeld calculation of the linear part of the self-ﬁeld (σi,Fi

) shares

the similar expression as the second-order moments of spatial coordinates of the bunch

distribution [34, 44]. That is the linear approximation of the force distribution depends

63

largely on the second-order moments of the distribution, which is consistent with the

argument by Sacherer [38].

The Kapchinsky-Vladimirsky (K-V) envelope equations are well-established in the ac-

celerator community for describing both continuous and bunched beams. It was initially

proposed to model beams with uniform density proﬁle and later extended to continu-

ous beams with elliptical symmetry as well as bunched beams with ellipsoidal symme-

try. [38] The SSA model is mathematically equivalent to the K-V envelope equations as

the SSA uses similar rms (root-mean-square) bunch statistics for uniform bunch evolu-

tion as Sacherer mentioned in his paper. The only diﬀerence is that the K-V envelope

equation combines the time-derivatives of bunch statistics in each direction into one or-

dinary diﬀerential equation. Speciﬁcally from the following two of the bunch statistics

dynamics [44]:

γi

(cid:32)

dσ2
i
dt
dγi
dt

=

=

2
m
1
m

ηi +

γ2
i
σ2
i

(cid:33)

+ σi,Fi

dσi
dt

=

1
2σi

dσ2
i
dt

=

1
m

γi
σi

(3.61a)

(3.61b)

(3.62)

Since we also have

64

we can rewrite Eq. 3.61a as:

(cid:33)
(cid:33)

− γi
σ2
i

dσi
dt

ηi +

+

γ2
i
σ2
i

(cid:35)

σi,Fi
σi

− γi
σ2
i

1
m

γi
σi

(cid:32)
(cid:34)

d2σi
dt2 =

=

=

=

1
m

1
m

1
σi

dγi
dt

(cid:32)

1
mσi

ηi

m2σi
ε2
i

m2σ3
i

+

+

σi,Fi
mσi
σi,Fi
mσi

The K-V envelope equations can be reproduced by the SSA model as:

dt2 − ε2
d2σi

i

m2σ3
i

/m

− σi,Fi
σi

= 0

(3.63)

(3.64)

The SSA model counterpart related to the varying emittance in K-V envelope equations

will be discussed next.

3.4.4 SSA with emittance evolution

The SSA model is capable of estimating the bunch evolution with varying emittance. We

can rewrite the SSA model (Eq. 3.6c) as:

(cid:32)

(cid:33)

dηi
dt

=

d
dt

ε2
i
σ2
i

= −2γiηi
mσ2
i

+

dε2
i
dt

1
σ2
i

(3.65)

The ﬁrst term on the right-hand side (RHS) can be considered as the energy transfer from

the linear local average motion to the local momentum ﬂuctuation. The second term on

the RHS represents the energy transfer from the potential energy, describing the impact

65

of a varying emittance on local momentum ﬂuctuation. The detailed form of these two

terms can be found in Eq. 3.6. Currently, in our research, the emittance evolution can be

obtained numerically from simulation data. We discuss our preliminary thoughts about

the emittance evolution process in the next chapter. A more thorough study discussing

the corresponding analytical models will be presented in future research.

Furthermore, the uniform density proﬁle evolution is not perfectly self-similar due to

the presence of a non-zero emittance. Speciﬁcally, the edge of the uniform distribution

will develop a non-uniform tail because of the local momentum spread close to the edge.

Therefore, even for the uniform distribution, we might expect non-conserved emittance

evolution.

As we present all the models, we test their performance and limitation by comparing

them with N -particle simulations, especially the longitudinal focusing of the electron

bunches in the next chapter.

66

Chapter 4

Numerical Simulations

In this chapter, we present the molecular dyamics (MD) simulation results for the evo-

lution of uniform ellipsoidal electron bunches focused by linear longitudinal chirps. We

discuss the comparison between the SSA model predictions and the MD simulation results

for further understanding of the electron bunch evolution (especially around crossover)

and the limitations of the SSA model.

4.1 Simulation method

Although there are many simulation methods that run faster than MD in studying the

electron bunch evolution with Coulomb eﬀects, such as the Particle-In-Cell (PIC) method,

we rely on the MD simulation due to the fact that it fully preserves the discrete aspect

of the Coulomb interaction information between electrons. Our simulation starts from a

uniform ellipsoidal bunch with a non-zero initial emittance, which we refer to as a “warm”

initial condition. Then, we solve non-relativistic equations of motion for every electron

using the velocity-Verlet integration. The Coulomb interaction is solved using the Fast

Multipole Method (FMM) from the fmmlib3d library [45] to calculate the ﬁeld at the

spatial position of each electron.

We present a simple overview about our MD simulation before we discuss the warm

67

initial conditions and the simulation results. The simulation is carried out in a rest frame

to simulate the evolution of an electron bunch in its center of mass (COM) frame. Each

particle represents one individual electron. Therefore, each particle/electron j has six

coordinates in the COM frame to describe the phase space position at time t during

the simulation: (rj = xj, yj, zj) for its spatial position in the COM frame and (vj =

vxj, vyj, vzj) for its relative velocity with respect to the COM.

Second, the Coulomb interaction on particle j from all other electrons is calculated by

the Fast Multipole Method (FMM), a numerical technique that speeds up the calculation
of long-range forces such as Coulomb interaction [45]. It achieves O(N ) time complexity

(where N is the total number of particle/electron) by expanding the corresponding Green’s

function using a multipole expansion, which equivalently groups sources that lie close

together and treats them as if they are a single source for long-range interaction.1 This

is used to calculate the ﬁeld to a speciﬁed accuracy without the use of a spatial mesh.

Third, we use the velocity-Verlet integration scheme as our non-relativistic pusher

to propagate the bunch evolution forward in time. Velocity-Verlet integration, similar

to the leapfrog integration, is a second-order integration method, that updates positions

and velocities at interleaved time points, staggered in such a way that they “leapfrog”

over each other to enhance accuracy into second-order using only one additional ﬁeld

solve calculation. With the force information Fj(t) from the FMM calculation based on

the position of particles at time t, the corresponding non-relativistic equation of motion

1An overview of the algorithm and implementation of FMM can be in this lecture notes by Beatson

and Greengard [46].

68

(EOM) for each individual particle (j) is then:

(cid:18)

rj(t + ∆t) = rj(t) +

vj(t) +

(cid:19)

· ∆t

Fj(t)
me

∆t
2

vj(t + ∆t) = vj(t) +

Fj(t) + Fj(t + ∆t)

me

· ∆t
2

The detailed implementation scheme for this integration is then:

(4.1a)

(4.1b)

1. propagate the velocity to the middle-point: vj(t + 1

2∆t) = vj(t) +

Fj (t)
me

∆t
2 using

the force calculation from the prior step

2. propagate the position to the next time-point: rj(t + ∆t) = rj(t) + vj(t + ∆t

2 )∆t

3. Calculate the force Fj(t + ∆t) based on the new charge distribution r(t + ∆t)

4. propagate the velocity to the next time-point: vj(t+∆t) = vj(t+ ∆t

2 )+

Fj (t+∆t)

me

∆t
2

Essentially, the velocity-Verlet takes advantage of staggering the calculation time point

for position and velocity to apply mid-point estimation for both the position and the

velocity evolution at once.

To assure the validity of the simulation, we check the conservation of energy by calcu-

lating the kinetic energy and potential energy of the bunch regularly during the simulation.

In this chapter, we present the MD simulation data by one line for the mean value of the

90 sample runs and one same-colored area shows the size of standard deviation about the

mean across the samples. A detailed discussion about the 90 samples will be present in

the next section.

The simulations were conducted using in-house code with the ﬁeld calculation done

by the fmmlib3d library. This code has been validated through comparison to another in-

69

house code implementing the brute-force MD code accelerated by GPGPU parallelization

using CUDA [47].

4.2 Warm initial condition

In this section, we present the “warm” initial condition for the equation of motion (Eq.

4.1) in our MD simulations.

In the SSA model and MD simulation, we assume the

longitudinal compression process for the probing electron bunch takes place during the

free-drifting region between the RF cavity and the specimen as showed in Fig. 1.1. The

warm initial condition is a thermalized bunch of particles that represents an electron bunch

with a non-zero emittance as our attempt to mimic the bunch phase-space distribution

right after the RF cavity. The warm initial condition helps us to prepare a phase-space

distribution of an electron bunch with non-zero emittance.

We ﬁrst place about 19100 electrons inside a simulation box with periodic boundary

condition (PBC) at the target density, which is 10, 000 electrons for a prolate spheroid

with the semi-axes of (10µm, 10µm, 75µm). The starting position of each electron is

randomly drawn from a uniform distribution and the starting momentum is zero.

Then the electrons interact with each other until the bunch reaches equilibrium. We

set the thermalization time to be over 10 plasma oscillation periods because (1) we ﬁnd

that this time period is long enough for bunches to reach equilibrium, and (2) the proper

travel time from the photo-cathode to the RF cavity is close to the duration of the

thermalization process. We refer to this process as the sample thermalization. Since we

employ the periodic boundary conditions during the sample thermalization, the Particle-

Particle-Particle-Mesh (PPPM) method is the preferred standard method as the Particle-

70

Mesh part of PPPM takes advantage of the highly eﬃcient Fast Fourier Transform (FFT)

method [48] to speed up the evaluation of the long-range interaction from the PBC, while

the normal FMM is not as eﬃcient for this particular situation as PPPM. Therefore, the

thermalization process of the initial conditions are prepared with the high-performance

PPPM feature in LAMMPS [49].

Figure 4.1: The z − pz phase-space distribution of a warm initial condition. The initial
local momentum ﬂuctuation has a Gaussian proﬁle and the corresponding linear chirp is
added based on the spatial coordinates of each electron.

At the end of the thermalization process, we select the electrons which are inside the

desired prolate spheroidal region to construct one sample of the warm initial condition.

One sample of the warm initial conditions is shown in Fig. 4.1. As the initial position is

random, the resulting samples will have slight diﬀerences in the total number of particles.

71

To mitigate the randomness of starting position, we prepare 90 such samples to obtain

reasonable statistical analysis for simulation results. To apply diﬀerent initial linear chirp

into the warm initial conditions, the momentum of each particle is adjusted based on its

spatial coordinates as presented in Eq. 3.5.

4.3 Width evolution

In this section, we present the comparison between the MD simulation results and the SSA

model predictions for the bunch evolution, speciﬁcally, the longitudinal focusing process.
Several simulations were performed for initial linear chirps: −0.7ω∗
where ω∗

c is obtained from the MFT calculation in Sec. 2.1.4. As shown in Fig. 4.2, the

c , −1.0ω∗

c and −1.5ω∗
c ,

SSA model predictions deviate from the MD simulation results in three familiar aspects:

(1) a slightly larger minimum width, (2) the critical time to reach the minimum width is

smaller and (3) the discrepancy between SSA and MD is most signiﬁcant when the bunch

is compressed at the critical chirp, where the bunch evolution is most sensitive to diﬀerent

emittance as we discussed in the last section. We previously saw similar trends in the

minimum width and the time to reach the minimum width as we discussed in Fig. 3.2.

As the input for the SSA model prediction in Fig. 4.2 already includes the mean value of

the initial emittance across the 90 simulation samples, this suggests that the longitudinal

emittance may be decreasing during the MD simulations. The rms emittance evolution

in panel (b) of Fig. 4.4 conﬁrms the decrease of longitudinal emittance, with several

interesting features we are going to discuss later.

As discussed previously, the mathematical mechanism behind the conservation of emit-

tance in the SSA is the vanishing of the following term associated with emittance evolution

72

Figure 4.2: Comparison of longitudinal width evolution with diﬀerent initial chirps from
SSA and MD simulations. For each initial chirp, the connected-dots correspond to the
average of the 90 MD simulations, the thick solid line is the SSA model prediction with
the average emittance of the warm initial conditions. The thick dotted line is the modiﬁed
SSA considering the emittance evolution in Fig. 4.4. The discrepancy between SSA and
MD simulations is driven by the longitudinal emittance decrease, seen in the simulations,
which is conﬁrmed by the good agreement between the modiﬁed SSA and simulations.

in Eq. 3.6:

σ2
i σpi,Fi

− σi,pi

σi,Fi

= 0

Currently we are working on the analytical theory to predict the evolution of these terms

in either a continuum approximation of the SSA model or the discrete particle settings

like the MD simulations.

However, we can extract the evolution of emittance from the simulation data to better

73

0100200300400500600700t(ps)0.00.20.40.60.81.0z(t)/z(0)FMM(0.7c)SSA(0.7c)SSAe(0.7c)FMM(1.0c)SSA(1.0c)SSAe(1.0c)FMM(1.5c)SSA(1.5c)SSAe(1.5c)capture the accuracy of the model predictions. Speciﬁcally, as we discussed in Sec. 3.4.4,

we replace Eq. 3.6c by Eq. 3.65:

dηi
dt

=

d
dt

(cid:33)

(cid:32)

ε2
i
σ2
i

= −2γiηi
mσ2
i

+

dε2
i
dt

1
σ2
i

in the SSA model with the time derivative of emittance square (dε2

i /dt). The ﬁrst term

on the right-hand side (RHS) can be considered as the coupling between the linear local

average momentum and the local momentum ﬂuctuation. This term comes from the non-

interacting kinematics as expansion reduces ηi and compression increases ηi. The second

term represents the energy transfer between the potential energy to ηi, which can be

obtained from the simulation data.

The evolution of the longitudinal width and the kinetic energy associated with the

electron motion in the longtudinal direction using this modiﬁcation can be seen as the

dotted lines in Fig. 4.2 and Fig. 4.3, respectively. Excellent agreement between these

modiﬁed SSA model predictions and the MD simulation results suggets that varying

emittance is the main factor causing the discrepancy between the longitudinal spatial

variance and longitudinal kinetic energy evolution of the (constant-emittance) SSA model

and the MD simulations. Therefore, if the evolution of the covariance term σpi,Fi

and

σi,Fi

can be understood and modeled, we should be able to obtain the SSA model that

captures the expected behavior of electron bunches to a high degree of accuracy.

74

Figure 4.3: Comparison of kinetic energy associated with the longitudinal motion of the
electrons relative to the bunch center of mass with diﬀerent initial chirps from SSA and
MD simulations. For each initial chirp, we have the connected-dots are the corresponding
90 MD simulations, the thick solid line for SSA prediction with the average emittance of
those warm initial conditions and the thick dotted line for the modiﬁed SSA considering
the emittance evolution in Fig. 4.4.

4.4 Emittance evolution

In this section, we discuss our preliminary thoughts about the emittance evolution process,

speciﬁcally the electron bunch focusing process. More thorough studies and corresponding

analytical models will be presented in future studies.

75

For now, the emittance evolution is obtained numerically from the simulation data

as shown in Fig. 4.4. The emittance evolution in Fig. 4.4, especially the simultaneous

Figure 4.4: (rms) Emittance evolution in (a) transverse εx and (b) longitudinal εz direc-
tion for three diﬀerent initial chirps.

increase for both the longitudinal and transverse emittance, cannot be explained by a

heat transfer mechanism employed in the literature [50]. On one hand, as can be seen in

panel(b) of Fig. 4.4, the longitudinal emittance increases slightly at the beginning of the

simmulation followed by a gradual decrease. In addition, there is another rapid increase

in the longitudinal emittance close to tc for the simulations in the crossover regime. On

the other hand, in panel (a) of Fig. 4.4, the transverse emittance also has a slight increase

at the beginning of the simulation followed by a gradual increase. Notice that again for

76

0.00.51.01.52.0x(nm)Transverse(a)FMM(0.7c)FMM(1.0c)FMM(1.5c)0100200300400500600700t(ps)0.80.91.01.11.21.3z(nm)(b)Longitudinalthe simulations in the crossover regime, there is a rapid increase around t = tc as well.

We would like to emphasize that the rapid increase in the emittance of both directions are

almost coincidental — an observation that is not predicted in the literature and cannot

be explained by the heat transfer mechanism.

We propose two mechanisms, Disorder-Induced Heating (DIH) and emittance transfer,

for driving the emittance evolution in the MD simulations [33] as discussed in the following

two subsections.

4.4.1 Disorder-induced heating (DIH)

DIH in the plasma community describes the heating process during the transition of

a bunch from a disordered state to an ordered state which is structured by Coulomb

forces [51–55]. We argue that there are two phases during the MD simulations where

signiﬁcant emittance growth is generated by this mechanism: (1) the sudden removal of

the thermalization conﬁnement at the beginning of the simulation and (2) around the

critical time tc, especially for the simulations in the crossover regime.

First, we attribute the emittance increase at the beginning of the simulation to DIH.

With the sudden removal of surrounding electrons and the periodic boundary condition,

the already thermalized electron bunch is suddenly not at equilibrium in the new cir-

cumstance. The bunch can be considered as a disordered state and it starts to evolve

towards its new equilibrium with the inﬂuence of Coulomb interaction. The heat transfer

from potential energy to the local momentum ﬂuctuation during this initial phase of the

simulation can then be explained by DIH.

The second situation for DIH induced emittance growth happens around the critical

77

time (tc) or focal point where a bunch reaches its minimum. In Fig. 4.5, we present the

bunch density distribution in MD simulations throughout the crossover evolution, and

its comparison with the density proﬁle of a uniform ellipsoid with the same longitudinal

spread. Since the diﬀerence is most signiﬁcant at tc, DIH suggests a signiﬁcant amount of

heat is injected into the kinetic energy associated with the local momentum ﬂuctuations.

In other words, the focusing process of the uniform bunch (especially for the cases in

crossover regime) pushes the bunch further away from its equilibrium state, resulting

in a highly non-equilibrium state at tc. The relaxation of this highly non-equilibrium

state is another demonstration of disorder induced heating [33]. As such DIH coming

from the potential energy is isotropic, we see a sudden increase in the emittance in both

the longitudinal and transverse direction almost simultaneously. A detailed discussion

about this DIH can be found in the following sections. Moreover, Eq. 3.65 also indicates

this emittance increase is further ampliﬁed at crossover where σ2

z reaches its minimum.

Therefore, a signiﬁcant longitudinal emittance boost is observed around tc in the crossover

case in panel (b) of Fig. 4.4.

4.4.2 Emittance transfer between degrees of freedom

In our MD simulations for the longitudinal focusing of the electron bunch, the longitudinal

emittance decreases due to the heat transfer to the transverse direction. Here the expres-

sion 1
2me

ηi is considered to be the SSA temperature of each degree of freedom because (1)

1

2me

ηi can be considered as the kinetic energy or temperature associated with the local

momentum ﬂuctuation, and (2) ηi is the variance of the local momentum ﬂuctuations

with respect to the linear local average momentum.

78

Figure 4.5: Longitudinal density proﬁle evolution of crossover case (−1.5ω∗
c ) with the
longitudinal width scaled by its standard deviation. In each panel, solid line represents
the average of simulation samples and dotted line is the density proﬁle for an equivalent
uniform spheroid of the same longitudinal width as the simulation data. The deviation
from uniformity is most prominent at crossover. However, tails of the distribution are
present in simulations for all initial chirps.

79

0.000.050.100.150.200.250.300.35(a)FMM@0.50tcUniform0.000.050.100.150.200.250.300.35(b)FMM@0.75tcUniform0.00.10.20.30.4z(1/z)(c)FMM@1.00tcUniform0.00.10.20.30.4(d)FMM@1.25tcUniform3210123z (z)0.000.050.100.150.200.250.300.35(e)FMM@1.50tcUniformFigure 4.6: The variance of the local momentum ﬂuctuations for the case in bounce-back
regime (0.7ωc) on a Log-Log scale.

Although the equilibrium requirement for a strict deﬁnition of temperature is not

appropriate for the highly non-equilibrium electron bunches in this analysis, we want to

borrow the temperature concept to facilitate our discussion about the heat/emittance

transfer. At the begining of the simulation, the initial value of the SSA temperature in

both longitudinal and transverse directions are similar to each other. Eq. 3.6c and Eq.

3.65 suggest that the longitudinal SSA temperature (ηz) increases as the electron bunch

focuses in the longitudinal direction, while the transverse SSA temperature decreases due

to the expanding transverse direction. Therefore, this SSA temperature diﬀerence between

the two direction increases as the bunch evolves towards the focal point. We would expect

80

a heat transfer from the hotter longitudinal direction to the cooler transverse direction

in such circumstances. This heat transfer is achieved through the exchange between

potential energy and the variance of local momentum ﬂuctuation ηi of each direction.

The immediate result of this heat transfer can be seen in Fig. 4.6 where ηz from the MD

simulations is smaller than that from the SSA model while ηx in the simulation is larger

than the model prediction.

This heat transfer is further illustrated in Fig. 4.7 as we plot the term associated with

the impact of a non-conserved emittance in Eq. 3.65 for the dynamics of ηi. After the

initial DIH eﬀect at the beginning, the increasing SSA temperature diﬀerence generates

stronger heat transfer from longitudinal direction to transverse direction. Prior to tc,

Figure 4.7: The heat transfer from potential to ηi (a) and (b) for bounce-back case
(0.7ωc); (c) and (d) for cross-over case (1.5ωc). The red dotted lines are the time when
the transfer in longitudinal direction switches sign.

the consistent and accelerating heat transfer from the longitudinal direction into the

transverse direction is evidence of the emittance transfer between the two degrees of

81

freedom.

Around tc, a heat inﬂux from the potential energy is observed in both the longitudinal

and transverse directions as the electron bunch relaxes from the highly non-equilibrium

state around tc. Notice that the heat inﬂux for the bounce-back cases are smaller than

the cross-over cases, and the non-equilibrium eﬀect is more drastic for the cross-over cases

as the density proﬁles are further away from unifrom than the bounce-back case. The

discussion about the density proﬁle at tc can be found in next section. Furthermore, the

non-equilibrium state at tc for cross-over case is so drastic that the heat continues to

transfer from the potential energy to the local momentum ﬂuctuation well beyond tc, as

can be seen in panel (c) and (d) in Fig. 4.7.

4.4.3 Phenomenological description of the emittance evolution

Due to the similarity in the longitudinal width evolution of the MD simulation results and

the (constant-emittance) SSA predictions in Fig. 4.2, the emittance evolution is largely

driven by the change in ηi.

The two competing mechanisms result in the following three phases for the emittance

evolution for bunches that are compressed longitudinally: (1) Initial phase dominated by

DIH: at the beginning of simulation, the DIH due to the removal of the thermalization

conﬁnement is the driving factor that increases the emittance in both longitudinal and

transverse direction; (2) Emittance transfer phase: the SSA temperature diﬀerence in-

creases as the bunch evolves, the emittance transfer mechanism overwhelms the initial

DIH, resulting in the decrease in longitudinal emittance and increase in transverse emit-

tance; (3) DIH around tc:

for the cases in the cross-over regime, the bunch reaches a

82

highly non-equilibrium state when it approaches its minimum width around tc, creating

another circumstance for a signiﬁcant DIH eﬀect. The relaxation of this DIH leads to the

sudden emittance growth in both the longitudinal and transverse direction in Fig. 4.4.

4.5 Density proﬁle at tc

In this section, we further study the highly non-equilibrium state and the corresponding

DIH eﬀect around tc. We present the spatial density proﬁle of one MD simulation sample

in the cross-over regime and the corresponding SSA estimation and uniform density proﬁle

in Fig. 4.8.

We generate the histogram and the corresponding kernel density estimation (KDE)

from the z coordinates of the electrons in the bunch at tc. By comparing the simula-

tion result with the uniform density proﬁle of the same longitudinal variance, it is clear

that the density proﬁle in the longitudinal direction of the MD simulation at tc deviates

signiﬁcantly from the uniform proﬁle.

4.5.1 Density proﬁle prediction from SSA and its linear chirp

assumption

We discuss the density proﬁle at tc determined from the SSA model prediction, followed

by a discussion of its linear-chirp assumption.

The density proﬁle at any given time t > 0 can be obtained by progressing the phase-

space distribution from the known initial distribution of (˜z, ˜pz) before adding the chirp

83

Figure 4.8: Longitudinal density proﬁle at tc of one crossover case (−1.5ω∗
c ). The grey
histogram represents the electron distribution along the longitudinal direction of the sim-
ulation sample at tc. The green solid line is the corresponding kernel density estimation
(KDE) of the density distribution. The blue solid line is the corresponding density pro-
ﬁle from the SSA calculation of that sample. The orange dotted line is representing the
density distribution of a uniform spheroid with the same longitudinal spread.

to (z(0), pz(0)):

dQ = S(˜z)d˜zG( ˜pz)d ˜pz

(4.2)

with S(˜z) representing a uniform distribution with the spheroidal shape and G( ˜pz) repre-

senting an independent Gaussian proﬁle of initial momentum spread without a chirp, as

84

shown in Fig. 4.1:

S(˜z) = πa2

G( ˜pz) = exp

(cid:34)
(cid:20)

(cid:18) ˜z
(cid:21)

c

2

(cid:19)2(cid:35)

(4.3a)

(4.3b)

1 −

− ˜pz
2η

from the initial bunch size (a(0), c(0)). Adding the linear chirp γ/σ2 simply means one

extra transformation:

 z(0)

 =

 1

pz(0)

γ/σ2 1

 ·

 z0



0

pz0

The transfer matrix T (t) for the trajectory of every particle according to the linear force

in SSA is then:

 z(t)

pz(t)

 = lim

∆t→0



t/∆t(cid:89)

i=0

1

1
me

∆t

KF z[α(i∆t)] · ∆t

1

 ·

 z(0)

 = T (t) ·

 z(0)



pz(0)

pz(0)

(4.5)

where T (t) and its inverse T−1(t)can then be obtained through numerical solutions for

the dynamics of the SSA model:

 z(0)
 = T−1(t) ·
 z(t)
 =
T−1
 =
 1
 =
 z(0)
 ·
 1

pz(0)

pz(t)

0
−γ/σ2 1

pz(0)

0
−γ/σ2 1

11 z(t) + T−1
T−1
21 z(t) + T−1

12 pz(t)

22 pz(t)

 · T−1(t) ·


 z(t)

pz(t)



and

 ˜z

˜pz

(4.4)

(4.6)

(4.7)

85

So then the longitudinal density distribution function of the SSA model reads:

ρ(z; t)dz =

S(˜zz(t),pz(t))d˜zG(˜zz(t),pz(t), ˜pzz(t),pz(t))d ˜pz

(cid:16)

−

˜pzz(t),pz(t)

2ηz

(cid:17)2

 dpz

(4.8)

(cid:90)

(cid:90)

(cid:34)

(cid:18) ˜zz(t),pz(t)

(cid:19)2(cid:35)

c

exp

= dz

πa2

1 −

with the help of the Liouville’s theorem, i.e. d˜zd ˜pz = dz(t)dpz(t). Therefore, we obtain

the projected density proﬁle according to the SSA model as the blue solid line in Fig. 4.8,

which also deviates quite signiﬁcantly from a uniform density proﬁle as the SSA model

assumed.

This derivation has a caveat that the transformation leads to a bell-shaped density

proﬁle while the driving T matrix is calculated following the uniform density proﬁle as-

sumption in the SSA model.

In other words, we model the bunch evolution using a

uniform density proﬁle while the actual density distribution is not strictly uniform. How-

ever, from the accuracy prospective, if the deviation from the uniform density proﬁle is

not too signiﬁcant (for bounce-back case) and not lasting too long (for cross-over case),

the predictions of SSA are still in a good agreement with the simulation results as we

can see from Fig. 4.2. The cases driven by the critical chirp present the most signiﬁcant

discrepancy between the simulation and the model.

This deviation from uniform density proﬁle illustrated both the limitation and versa-

tility of the SSA model. The limitation lies as it relies on the linear chirp assumption to

conserve emittance while only the uniform bunch with zero emittance can maintain that

linear chirp assumption throughout the bunch evolution. The versatility of SSA is that

with the addition of the emittance change term, the modiﬁed SSA model is a complete

86

statistical model that focuses on the evolution of bunch statistics.

4.6 Bunch phase space display for cylindrical bunch

The MD simulation provides the detailed phase-space structure for the bunch evolution.

It can be particularly helpful for the bunch proﬁles that are hard to model. For example,

the initial bunch shape might be closer to a cylinder after passing through an aperture,

which is helpful in reducing the emittance by removing the hot region on the out-skirts

of the bunch [17]. For modeling the focusing process of such bunches, we simulate a uni-

form bunch starting with cylinder shape with a linear chirp in the longitudinal direction.

Fig. 4.9 shows an interesting non-linear phase-space and real-space electron distribution.

Figure 4.9: The longitudinal phase-space (ﬁrst row) and real-space (second row) distri-
bution evolution of a uniform bunch, which starts with a cylindrical shape, zero emittance
and a linear longitudinal chirp. Electrons are colored based on their initial longitudinal
position.

Speciﬁcally, by comparing the phase-space and real-space distribution, we realized that

87

the non-linearity in the phase-space distribution (the hook on the both ends) comes from

the center of the front and back end in the longitudinal direction, as the net eﬀect of

Coulomb repulsion from the “extra” electrons that are initially outside of the ellipsoidal

region is most signiﬁcant in those two regions. Furthermore, the Coulomb repulsion from

those extra electrons prevents the electrons in those two regions from crossing the bunch

center, resulting in a partial-crossover situation.

4.7 Discussion

In this chapter, we present the MD simulation results for the longitudinal focusing pro-

cess of uniform spheroidal electron bunches. We equilibrate the electrons for the warm

initial condition to minimize the initial emittance growth due to DIH. The comparison

of the longitudinal width evolution between the MD simulations and the SSA predictions

shows the impact of a varying emittance on bunch evolution. We propose two compet-

ing mechanisms for the change of emittance throughout the compression process. The

DIH increases the emittance in both degrees of freedom while the diﬀerence in the SSA

temperature generates emittance transfer between degrees of freedom. In addition, the

non-uniform density proﬁle at the focal point (t = tc) introduces signiﬁcant DIH that in-

creases the emittance in both longitudinal and transverse directions. The MD simulations

also unveil other interesting phase-space structures during the bunch focusing process.

In this section, we extend our discussion about the MD simulation results and the

comparison between the MD simulation and the SSA model.

88

4.7.1 More on the emittance evolution

In Sec. 4.4.2, we discussed the heat transfer due to the SSA temperature diﬀerence in the

longitudinal and transverse directions. In addition to the bounce-back example in Fig.

4.6, we present the evolution of ηi for a crossover case, with the comparison to the SSA

model prediciton in Fig. 4.10. The heat inﬂux in the transverse direction is much more

Figure 4.10: The variance of the local momentum ﬂuctuations for the cases in the cross-
over regime (1.5ωc) on a Log-Log scale.

signiﬁcant than that in the bounce-back case. As the bunch evolves far beyond tc, we

would expect ηz and ηx to become closer due to the emittance transfer. However, that

is not obvious from Fig. 4.6 and Fig. 4.10. The possible explanation could be that the
exchange between ηi and the linear average motion within each direction, − 2γiηi
meσ2
i

, can

be signiﬁcantly larger than the heat transfer between directions, as the γi increases along

89

with the evolution.

While the two mechanisms for emittance dynamics, DIH and emittance transfer, is

qualitatively reasonable, it would be more convincing if we can provide more quantitative

estimation to diﬀerentiate the eﬀect of these two mechanisms. Speciﬁcally, our future

work will focus on analytical models to quantify the DIH for the relaxation process from

the bell-shaped density proﬁle to a uniform distribution as shown in Fig. 4.5. One of the

candidate is the Wrangler’s theorem [56–61], which estimates the emittance growth due

to the bunch density distribution evolving from a “less” uniform proﬁle towards a “more”

uniform one. Nonetheless, using these two mechanisms, we reach a phenomenological

description about the emittance evolution during the longitudinal focusing process of

uniform ellipsoidal bunches.

4.7.2 Sampling error

As we presented in Fig. 4.4, the emittance of a uniform bunch evolves quite dramatically

in our MD simulations. Another perspective of understanding the origin of emittance evo-

lution might be the sampling error. Speciﬁcally, the SSA model (and all other analytical

models) treats the electron bunch as a continuum object, while the electron bunch in MD

simulations consists of discrete electrons randomly drawn from the uniform distribution.

Thus, the SSA model can be considered as a statistical description for the population and

each simulation is just one ﬁnite-sized sample drawn from the population. In this sense,

the simulation introduces sampling errors so that the measurement of the sample deviates

from the true value of the population. For example, we assume that some bunch statistics

like (cid:104)z(cid:105), σz,δz are zero in the model, which are most likely non-zero in the simulations.

90

It is unclear if there is a recovery mechanism for this sampling error. The recovery

mechanism here is designated to the intrinsic feature of the statistics or kinematics or

Coulomb interaction that reduces the deviation between the sample and the population

towards zero as the sample/bunch evolves. From the model and existing simulation results,

there is no such eﬀect or the magnitude is not signiﬁcant enough, i.e. the deviation between

a single simulation and the SSA model does not decrease from the initial sampling error.

Further analytical study, such as the impact due to diﬀerent size of sampling error, is

needed to help us understand the origin of this emittance change mechanism.

4.7.3 Comparing with the experiment values

To put our MD simulation setup into perspective, here are the typical values of some key

operating parameters in the actual UED experiments [18]:

• Number of electrons in a single electron bunch (Ne): 106 — 1.7 × 107

• Beam velocity (ve): 164.35 µm/ps

• Relativistic Lorentz factor (γCOM ): 1.1957

• RF cavity to specimen distance(L): 0.425m

• Pre-RF width (σ0) for Ne = 106 bunches: 5.8ps

• Pre-RF chirp (az) for Ne = 106 bunches: 0.00041ps−1

In the UED experiments, the RF focusing power is tuned to minimize the longitudinal

spread of the bunch when it arrives the specimen location.

91

For the bunch with Ne = 106, the number of electrons reduces to 104 passing an aper-

ture with diameter of 100µm after the RF cavity2. The corresponding bunch parameters

in the COM frame then reads:

• transverse radius: 50µm

• longitudinal variance in COM: σ2

z

(cid:12)(cid:12)t=0 =

(cid:17)2

(cid:16) σ0·ve

γCOM

= (797.2µm)2

• time from RF cavity to specimen in COM frame: τCOM = L/ve
γCOM

= 2.1ns

Our simulations are studying a density about 250 times higher than this set of experi-

mental parameters and a slightly diﬀerent initial aspect ratio. However, similar to our

argument in Sec. 2, the space-charge eﬀect we observed in the MD simulations is similar

for diﬀerent density, with the diﬀerence captured by a time-scaling factor Ω0. Therefore,

in order to compare with the experiments, we need to adjust the evolution of the simula-

tion by the corresponding Ω0. Judging from the comparison between the corresponding

SSA evolution estimation [44], the applied focusing chirp in the experiments ranges across

both the bounce-back and cross-over regimes.

2Discussion with Faran Zhou

92

Chapter 5

Conclusion

As mentioned in the introduction, this dissertation is centered on the analytical and nu-

merical study of the ellipsoidal electron bunch dynamics. We have examined the crossover

of electron bunches with uniform ellipsoidal proﬁles focused by a linear chirp as is typical

of the propagation of a probing electron bunch in ultrafast electron diﬀraction/microscope

systems.

We employed several analytic models to study the space charge eﬀects on the bunch

dynamics, the ﬁrst of which is an extension of the mean-ﬁeld model with ordinary dif-

ferential equations. Analysis of this mean-ﬁeld model leads to the identiﬁcation of a

longitudinal critical chirp. This critical chirp separates two regimes for particle trajecto-

ries in this model: bounce-back, where the particles reverse their direction at the waist of

the focusing process, and cross-over, where the bunch experiences a singularity where the

bunch width reduces to zero. We showed that time can be scaled by the initial plasma

frequency, and the critical chirp becoming dimensionless and solely depend on the initial

aspect ratio. The evolution of bunches with the same initial aspect ratio diﬀer only by

the time scale determined by the bunch’s plasma frequency [33]. The major drawback

of this model is that it requires a zero-emittance bunch, which limits the applicability of

this model comparing to the statistical methods.

We also presented the SSA model, a statistical method describing the second order

93

moments to model the evolution of an ellipsoidal electron bunch and explained how it

can reproduce the MFT by taking a zero emittance. We also discussed the linear chirp

assumption, explaining how it is the key assumption to the emittance conservation during

the bunch evolution according to the SSA model. In addition, we also shed light on the

statistical nature of emittance, noting that the force impact by the uniform density proﬁle

and the Gaussian proﬁle are close to each other. The expression for the force impact from

a changing emittance is captured by the additional term in the modiﬁed SSA model.

The modiﬁed SSA model is then equivalent to the K-V envelope equation. Through

the derivation of the SSA model, we pointed out that the application of the statistical

methods can extend beyond the uniform ellipsoidal bunch, while the accuracy of the SSA

prediction is mainly related to the discrepancy between the actual density proﬁle and the

uniform density proﬁle.

We then presented the MD simulation results for the longitudinal focusing process of

uniform spheroidal electron bunches. We equilibrate the electron bunches for the warm

initial conditions to minimize the initial emittance growth due to DIH. The comparison

of the longitudinal width evolution between the MD simulations and the SSA predictions

shows the impact of a varying emittance on bunch evolution. We propose two compet-

ing mechanisms for the change of emittance throughout the compression process. The

DIH increases the emittance in both degrees of freedom while the diﬀerence in the SSA

temperature generates emittance transfer between degrees of freedom. In addition, the

non-uniform density proﬁle at the focal point (tc) introduces signiﬁcant emittance growth

in both the longitudinal and transverse directions.

Moving forward, the next steps in the development of analytical modeling of the

94

probing electron evolution will focus on the theoretical understanding of the emittance

evolution, especially the quantitative estimation of the DIH and emittance transfer eﬀects

during the focusing process.

95

APPENDIX

96

Appendix

Detailed derivation for electrostatic

potential for uniform ellipsoid

In this appendix, we present the detailed derivation for Coulomb potential for uniformly

charged ellipsoids. The derivation is generally followed MacMillan’s book [39].

The surface of a given uniform ellipsoid with semi-axes (a, b, c) is deﬁned by equation:

ξ2
a2 +

η2
b2 +

ζ2
c2 = 1

(.1)

Let the interior point for which the potential is to be computed be P (x, y, z). On taking

P as the origin of a spherical coordinates system ρ, ϕ, θ with the transformation:

ξ = x + ρ cos ϕ cos θ

η = y + ρ cos ϕ sin θ

ζ = z + ρ sin ϕ

and the corresponding charge element

dq = ρcρ2 cos ϕdϕdθdρ

97

(.2a)

(.2b)

(.2c)

(.3)

with the volume charge density ρc = n · e (typically n is the electron number density).

Then, the electrostatic potential at point P (with respect to the potential zero point at

inﬁnity) is expressed as:

V (x, y, z) =

(cid:90)

E

ρc
4πε0

dq
ρ

= κ

(cid:90) + π

2

− π
2

(cid:90) 2π

(cid:90) ρ1(θ,ϕ)

0

0

We substitute κ = ρc
4πε0

to simplify the derivation.

ρ cos ϕdϕdθdρ

(.4)

The upper limit ρ1(θ, ϕ) of the integration with respect to ρ is a function of θ and ϕ,

since the integration is from P to a point on the surface of the ellipsoid. So, we can insert

Eq. .2a into Eq. .1 for ρ1 as:

Aρ2

1 + 2Bρ1 + C = 0

where

A =

B =

C =

cos2 ϕ cos2 θ

a2

x cos ϕ cos θ

a2
y2
b2 +

x2
a2 +

+

+

cos2 ϕ sin2 θ

b2

+

sin2 ϕ

c2

y cos ϕ sin θ

b2

+

z sin ϕ

c2

z2
c2 − 1

(.5)

(.6a)

(.6b)

(.6c)

With A being positive and C being negative for interior points, ρ1 is then the positive

root of Eq. .5 as:

−B +

ρ1 =

√
B2 − AC
A

(.7)

98

To compute V , we can take the advantage of(cid:82) ρ1

(cid:90) + π

2

(cid:90) 2π

− π
2

0

V =

κ
2

2B2 − AC − 2B
A2

2 ρ2

0 ρdρ = 1
√

B2 − AC

1, so that Eq. .4 becomes

cos ϕdϕdθ

(.8)

The second part of the integral vanishes as

(cid:90) + π

2

(cid:90) 2π

− π
2

0

R = κ

√

B

B2 − AC
A2

cos ϕdϕdθ = 0

(.9)

can be evaluated by taken in pairs the two elements canceling out each other, e.g. (θ0, ϕ0)
and (θ0 + π,−ϕ0). The integral reduces to
(cid:90) 2π

(cid:90) + π

2B2 − AC

2

cos ϕdϕdθ

(.10)

V =

κ
2

− π
2

0

A2

· x2
a2 +

cos2 ϕ sin2 θ

b2

· y2
b2 +

sin2 ϕ

c2

· z2
c2

zx cos ϕ sin ϕ cos θ

(cid:27) cos ϕdϕdθ

A2

(cid:27) cos ϕdϕdθ

c2a2

A2

(.11)

Substituting the B value from Eq. .6 and we get

V = κ [

2

(cid:26)cos2 ϕ cos2 θ
(cid:90) 2π
(cid:90) + π
(cid:26)xy cos2 ϕ cos θ sin θ
(cid:90) 2π
(cid:90) + π
(cid:90) + π
(cid:90) 2π

− π
2

− π
2

a2b2

a2

cos ϕdϕdθ

2

0

0

2

− π
2

0

A

]

+ 2

− C
2

+

yz cos ϕ sin ϕ sin θ

b2c2

+

99

By properly pairing the elements in the second integral vanishes, as

2

− π
2

(cid:90) + π
(cid:90) 2π
(cid:90) 2π
(cid:90) + π
(cid:90) 2π
(cid:90) + π

− π
2

0

2

0

2

− π
2

0

xy

yz

zx

cos2 ϕ cos θ sin θ

a2b2

cos ϕ sin ϕ sin θ

b2c2

cos ϕ sin ϕ cos θ

a2b2

· cos ϕdϕdθ

A2

· cos ϕdϕdθ

A2

· cos ϕdϕdθ

A2

with {(θ0, ϕ0), (2π − θ0, ϕ0)}

with {(θ0, ϕ0), (2π − θ0, ϕ0)}

with {(θ0, ϕ0), (π + θ0, ϕ0)}

we can have

(cid:18)cos2 ϕ cos2 θ
(cid:90) 2π
(cid:90) + π
(cid:90) 2π
(cid:90) + π

− π
2

cos ϕdϕdθ

a2

2

2

0

− π
2

0

A

V = κ

− κ
2

C

· x2
a2 +

cos2 ϕ sin2 θ

b2

· y2
b2 +

sin2 ϕ

c2

· z2
c2

(cid:19) cos ϕdϕdθ

A2

To further simplify the above expression, one may use the following:

(cid:90) + π

2

(cid:90) 2π

− π
2

0

W =

κ
2

cos ϕdϕdθ

A

(.12)

(.13)

The simpliﬁed form of the potential then reduces to:

V =

=

1
a

(cid:18) 1

∂W
∂a
∂W
∂a

a

x2 +

1
b
− W
a2

∂W
∂b

(cid:19)

y2 +

x2 +

1
c

(cid:18)1

b

∂W
∂c
∂W
∂b

(cid:19)

z2 − CW
− W
b2

y2 +

(cid:18)1

c

∂W
∂c

− W
c2

(cid:19)

(.14)

z2 + W

Since W is a function of the semi-axes (a, b, c), then so are all of its derivatives. Therefore,

the coeﬃcients of the quadratic term x2, y2 and z2 are functions of a, b, and c only. Also,

W is the potential at the center of the ellipsoid if taking P (x, y, z) = (0, 0, 0).

100

We can simplify the expression of W with the help of the following substitution:

M =

N =

cos2 ϕ
a2 +
cos2 ϕ
b2 +

sin2 ϕ

c2

sin2 ϕ

c2

(cid:90) 2π
(cid:90) π

0

2

dθ

M cos2 θ + N sin2 θ

sec2 θdθ

M + N tan2 θ

0

cos ϕdϕ

Eq. .13 then becomes

2

(cid:90) + π
(cid:90) π
(cid:90) π

− π
2

0

2

2

W =

κ
2

cos ϕdϕ

= 4κ

cos ϕdϕ

= 2πκ

0

= 2πκabc2

cos ϕdϕ√
M N

(cid:90) π

2

(cid:113)

0

(.15a)

(.15b)

(.16)

(.17)

(.18)

(a2 sin2 ϕ + c2 cos2 ϕ)(b2 sin2 ϕ + c2 cos2 ϕ)

To restore the symmetry within a, b, and c, we can introduce the following substitution:

sin ϕ =

c√
c2 + s

⇒ d sin ϕ
ds

= −

c

2(c2 + s)3/2

where s is the new variable of integration. Eventually, W can be written as:

(cid:90) ∞

0

W = πκabc

(cid:112)(a2 + s)(b2 + s)(c2 + s)

ds

101

with the term related to the derivative of W with respect to a reads:

1
a

∂W
∂a

− W

a2 =

πκabc

a2

= −πκabc

 ∂a(a2 + s)
(cid:90) ∞
(cid:90) ∞

∂a

− 1
2


(cid:112)(a2 + s)(b2 + s)(c2 + s)

1√
a2 + s

a2 + s

0

−

ds

1

0

(cid:112)(b2 + s)(c2 + s)

ds

(.19)

Therefore, the potential inside of a uniform ellipsoidal electron bunch may be expressed

as:

V (x, y, z) =

ρc
4ε0

abc

(cid:90) ∞

0

(cid:18)
1 − x2
a2 + s

(cid:19)

(cid:112)(a2 + s)(b2 + s)(c2 + s)

ds

− y2
b2 + s

− z2
c2 + s

(.20)

with zero point of potential at inﬁnity.

The corresponding electrostatic ﬁeld at the interior point P (x, y, z) is,

(cid:126)E(x, y, z) = Ex ˆx + Ey ˆy + Ez ˆz

(.21)

with ˆx, ˆy, and ˆz representing the unit vectors, respectively. The components of the ﬁeld

can be represented by the following equations:

Ex(x, y, z) = x · ρcabc
2ε0
Ey(x, y, z) = y · ρcabc
2ε0
Ez(x, y, z) = z · ρcabc
2ε0

(cid:90) ∞
(cid:90) ∞
(cid:90) ∞

0

0

0

ds

(a2 + s)(cid:112)(a2 + s)(b2 + s)(c2 + s)
(b2 + s)(cid:112)(a2 + s)(b2 + s)(c2 + s)
(c2 + s)(cid:112)(a2 + s)(b2 + s)(c2 + s)

ds

ds

(.22a)

(.22b)

(.22c)

102

BIBLIOGRAPHY

103

BIBLIOGRAPHY

[1] J. M. E. Portman, Ultrafast Science: A Multiscale Modeling Approach to Femtosecond
Electron Diﬀraction and Its Applications. PhD thesis, Michigan State University,
2014.

[2] 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, p. 111101, June 2005.

[3] Z. Tao, Development of High-Brightness Ultrafast Electron Microscope for Study-
ing Nanoscale Dynamics Associated with Strongly Correlated Materials. PhD thesis,
Michigan State University, 2014.

[4] R. D. Miller, “Mapping Atomic Motions with Ultrabright Electrons: The Chemists’
Gedanken Experiment Enters the Lab Frame,” Annu. Rev. Phys. Chem., vol. 65,
pp. 583–604, Apr. 2014.

[5] R. J. D. Miller, “Femtosecond Crystallography with Ultrabright Electrons and X-
rays: Capturing Chemistry in Action,” Science, vol. 343, pp. 1108–1116, Mar. 2014.

[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 tran-
sitions,” Proceedings of the National Academy of Sciences, vol. 101, pp. 1123–1128,
Feb. 2004.

[7] P. Baum, D.-S. Yang, and A. H. Zewail, “4D Visualization of Transitional Structures
in Phase Transformations by Electron Diﬀraction,” Science, vol. 318, pp. 788–792,
Nov. 2007.

[8] P. Baum and A. Zewail, “Femtosecond diﬀraction with chirped electron pulses,”

Chem. Phys. Lett., vol. 462, pp. 14–17, Sept. 2008.

[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 Diﬀraction,” Phys. Rev.
Lett., vol. 100, p. 155504, Apr. 2008.

[10] A. H. Zewail, “Four-Dimensional Electron Microscopy,” Science, vol. 328, pp. 187–

193, Apr. 2010.

104

[11] A. Feist, Next-Generation Ultrafast Transmission Electron Microscopy – Develop-
ment and Applications. PhD thesis, Georg-August-Universit¨at G¨ottingen, G¨ottingen,
Germany, 2018.

[12] 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 Nanocrys-
tals Revealed by Ultrafast Electron Crystallography,” Phys. Rev. Lett., vol. 104,
p. 123401, Mar. 2010.

[13] Z. Tao, H. Zhang, P. M. Duxbury, M. Berz, and C.-Y. Ruan, “Space charge eﬀects
in ultrafast electron diﬀraction and imaging,” Journal of Applied Physics, vol. 111,
p. 044316, Feb. 2012.

[14] 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.

[15] C. Kiseok, Development of a High-Brightness Electron Beam System towards Fem-
tosecond Microdiﬀraction and Imaging and Its Applications. PhD thesis, Michigan
State University, 2014.

[16] J. Portman, H. Zhang, K. Makino, C. Y. Ruan, M. Berz, and P. M. Duxbury, “Un-
tangling the contributions of image charge and laser proﬁle for optimal photoemission
of high-brightness electron beams,” J. Appl. Phys., vol. 116, p. 174302, Nov. 2014.

[17] B. S. Zerbe, X. Xiang, C.-Y. Ruan, S. M. Lund, and P. M. Duxbury, “Dynamical
bunching and density peaks in expanding Coulomb clouds,” Phys. Rev. Accel. Beams,
vol. 21, p. 064201, June 2018.

[18] J. Williams, F. Zhou, T. Sun, Z. Tao, K. Chang, K. Makino, M. Berz, P. M. Duxbury,
and C.-Y. Ruan, “Active control of bright electron beams with RF optics for fem-
tosecond microscopy,” Struct. Dyn., vol. 4, p. 044035, July 2017.

[19] A. A. Ischenko, I. V. Kochikov, and R. J. D. Miller, “The eﬀect of Coulomb repulsion
on the space-time resolution limits for ultrafast electron diﬀraction,” J. Chem. Phys.,
vol. 150, p. 054201, Feb. 2019.

[20] F. Zhou, Non-Equilibrium Phase Transformations in Charge-Density Wave and
Strongly Correlated Systems Studied by Coherent Femtosecond Electron Diﬀraction.
PhD thesis, Michigan State University, 2019.

105

[21] O. J. Luiten, S. B. van der Geer, M. J. de Loos, F. B. Kiewiet, and M. J. van der
Wiel, “How to Realize Uniform Three-Dimensional Ellipsoidal Electron Bunches,”
Phys. Rev. Lett., vol. 93, Aug. 2004.

[22] M. Grech, R. Nuter, A. Mikaberidze, P. Di Cintio, L. Gremillet, E. Lefebvre,
U. Saalmann, J. M. Rost, and S. Skupin, “Coulomb explosion of uniformly charged
spheroids,” Phys. Rev. E, vol. 84, Nov. 2011.

[23] J. Rosenzweig, A. Cook, R. England, M. Dunning, S. Anderson, and M. Ferrario,
“Emittance compensation with dynamically optimized photoelectron beam proﬁles,”
Nucl. Instrum. Methods Phys. Res. Sect. Accel. Spectrometers Detect. Assoc. Equip.,
vol. 557, pp. 87–93, Feb. 2006.

[24] B. J. Claessens, S. B. van der Geer, G. Taban, E. J. D. Vredenbregt, and O. J. Luiten,

“Ultracold Electron Source,” Phys. Rev. Lett., vol. 95, p. 164801, Oct. 2005.

[25] Y. Otake and B. J. Siwick, “Advanced diagnosis of the temporal characteristics of
ultra-short electron beams,” Nucl. Instrum. Methods Phys. Res. Sect. Accel. Spec-
trometers Detect. Assoc. Equip., vol. 637, no. 1, pp. S7–S11, 2011.

[26] S. B. van der Geer, M. J. de Loos, T. van Oudheusden, W. P. E. M. op ’t Root,
M. J. van der Wiel, and O. J. Luiten, “Longitudinal phase-space manipulation of
ellipsoidal electron bunches in realistic ﬁelds,” Phys. Rev. Spec. Top. - Accel. Beams,
vol. 9, p. 044203, Apr. 2006.

[27] P. Musumeci, J. T. Moody, R. J. England, J. B. Rosenzweig, and T. Tran, “Exper-
imental Generation and Characterization of Uniformly Filled Ellipsoidal Electron-
Beam Distributions,” Phys. Rev. Lett., vol. 100, p. 244801, June 2008.

[28] B. J. Siwick, J. R. Dwyer, R. E. Jordan, and R. J. D. Miller, “Ultrafast electron optics:
Propagation dynamics of femtosecond electron packets,” J. Appl. Phys., vol. 92, no. 3,
p. 1643, 2002.

[29] M. Kuwahara, Y. Nambo, K. Aoki, K. Sameshima, X. Jin, T. Ujihara, H. Asano,
K. Saitoh, Y. Takeda, and N. Tanaka, “The Boersch eﬀect in a picosecond pulsed elec-
tron beam emitted from a semiconductor photocathode,” Appl. Phys. Lett., vol. 109,
p. 013108, July 2016.

[30] Y. Zou, Y. Cui, M. Reiser, and P. G. O’Shea, “Observation of the Anomalous Increase
of the Longitudinal Energy Spread in a Space-Charge-Dominated Electron Beam,”
Phys. Rev. Lett., vol. 94, Apr. 2005.

[31] C. C. Lin, L. Mestel, and F. H. Shu, “The Gravitational Collapse of a Uniform

Spheroid.,” The Astrophysical Journal, vol. 142, p. 1431, Nov. 1965.

106

[32] P. Di Cintio, Dynamics of Heterogeneous Clusters under Intense Laser Fields. PhD

thesis, Technische Universit¨at Dresden, Dresden, Germany, 2014.

[33] X. Xiang, B. Zerbe, and P. M. Duxbury, “Crossover and the dynamics of uniform
electron ellipsoids focused by linear chirp,” working paper, Michigan State University,
2020.

[34] A. M. Michalik and J. E. Sipe, “Analytic model of electron pulse propagation in
ultrafast electron diﬀraction experiments,” J. Appl. Phys., vol. 99, no. 5, p. 054908,
2006.

[35] A. M. Michalik, E. Y. Sherman, and J. E. Sipe, “Theory of ultrafast electron diﬀrac-
tion: The role of the electron bunch properties,” J. Appl. Phys., vol. 104, no. 5,
p. 054905, 2008.

[36] 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, p. 124905,
Dec. 2010.

[37] Kapchinskij and Vladimirskij, “Limitations of Proton Beam Current in a Strong Fo-
cusing Linear Accelerator Associated With the Beam Space Charge,” in Proceedings,
(Geneva, Switzerland), pp. 274–287, CERN, 1959.

[38] F. J. Sacherer, “RMS Envelope Equations with Space Charge,” IEEE Trans. Nucl.

Sci., vol. 18, no. 3, pp. 1105–1107, 1971.

[39] W. MacMillan, The Theory of the Potential. MacMillan’s Theoretical Mechanics,

New York, NY: Dover, 1958.

[40] R. Gluckstren, “Scalar Potential for Charge Distributions With Ellipsoidal Symme-
try,” tech. rep., Physics Department, University of Maryland, Fermilab, May 1986.

[41] O. D. Kellogg, Foundations of Potential Theory. Berlin, Heidelberg: Springer Berlin

Heidelberg, 1967.

[42] Landau and Lifshitz, Electrodynamics Of Continuous Media, vol. 8 of Course of

Theoretical Physics. Pergamon Press, 1960.

[43] J. Buon, “Beam Phase Space and emittance,” tech. rep., Paris-11 Univ., 1992.

[44] B. S. Zerbe, X. Xiang, J. Williams, C.-Y. Ruan, and P. M. Duxbury, “Mean-ﬁeld
space charge eﬀects on longitudinally focussed electron bunches,” working paper,
Michigan State University, Nov. 2019.

107

[45] Z. Gimbutas and L. Greengard, “Computational Software: Simple FMM Libraries
for Electrostatics, Slow Viscous Flow, and Frequency-Domain Wave Propagation,”
Commun. Comput. Phys., vol. 18, pp. 516–528, Aug. 2015.

[46] R. Beatson and L. Greengard, “A short course on fast multipole methods,” Wavelets

Multilevel Methods Elliptic PDEs, vol. 1, pp. 1–37, 1997.

[47] R. G. Belleman, J. B´edorf, and S. F. Portegies Zwart, “High performance direct
gravitational N-body simulations on graphics processing units II: An implementation
in CUDA,” New Astronomy, vol. 13, pp. 103–112, Feb. 2008.

[48] B. A. Luty, M. E. Davis, I. G. Tironi, and W. F. Van Gunsteren, “A Comparison
of Particle-Particle, Particle-Mesh and Ewald Methods for Calculating Electrostatic
Interactions in Periodic Molecular Systems,” Molecular Simulation, vol. 14, pp. 11–
20, Dec. 1994.

[49] S. Plimpton, “Fast Parallel Algorithms for Short–Range Molecular Dynamics,” J.

Comput. Phys., vol. 117, p. 42, Mar. 1995.

[50] M. Reiser, Theory and Design of Charged Particle Beams. Weinheim: Wiley-VCH,

second ed., 2008.

[51] D. Gericke and M. Murillo, “Disorder-induced heating of ultracold plasmas,” Contrib.

Plasma Phys., vol. 43, pp. 298–301, Oct. 2003.

[52] M. S. Murillo, “Ultrafast Dynamics of Strongly Coupled Plasmas,” Phys. Rev. Lett.,

vol. 96, p. 165001, Apr. 2006.

[53] M. S. Murillo, “Ultrafast dynamics of neutral, ultracold plasmas,” Physics of Plas-

mas, vol. 14, p. 055702, May 2007.

[54] M. Lyon and S. D. Bergeson, “The inﬂuence of electron screening on disorder-induced

heating,” J. Phys. B: At. Mol. Opt. Phys., vol. 44, p. 184014, Sept. 2011.

[55] D. Murphy, R. E. Scholten, and B. M. Sparkes, “Increasing the Brightness of Cold
Ion Beams by Suppressing Disorder-Induced Heating with Rydberg Blockade,” Phys.
Rev. Lett., vol. 115, p. 214802, Nov. 2015.

[56] T. P. Wangler, K. R. Crandall, R. S. Mills, and M. Reiser, “Relation between Field
Energy and RMS Emittance in Intense Particle Beams,” IEEE Trans. Nucl. Sci.,
vol. 32, no. 5, pp. 2196–2200, 1985.

[57] J. Struckmeier and I. Hofmann, “The Problem of Self-Consistent Particle Phase Space

Distributions for Periodic Focusing Channels,” Part. Accel., vol. 39, p. 31, 1992.

108

[58] I. Hofmann and J. Struckmeier, “Generalized Three-Dimensional Equations for the
Emittance and Field Energy of High-Current Beams in Periodic Focusing Structures,”
Part. Accel., vol. 21, p. 30, 1987.

[59] O. A. Anderson, “Internal Dynamics and Emitiance Growth in Space Charge Domi-

nated Beams,” Part. Accel., p. 30, 1987.

[60] S. M. Lund and B. Bukh, “Stability properties of the transverse envelope equa-
tions describing intense ion beam transport,” Phys. Rev. ST Accel. Beams, vol. 7,
p. 024801, Feb. 2004.

[61] M. Reiser, “Free energy and emittance growth in nonstationary charged particle

beams,” Journal of Applied Physics, vol. 70, pp. 1919–1923, Aug. 1991.

109