ELFIN

Three Dimensional Free Electron Laser Simulation Code V.77

Takeshi Nakamura

SPring-8

Kamigori-cho, Ako-gun, Hyogo 678-12 JAPAN nakamura@spring8.or.jp

February 5, 1999

Contents

1 Abstract 2 2 INTRODUCTION 2 3 ASSUMPTION 2 4 FEL AMPLIFIER 4

4.1 FEL EQUATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 4.2 CALCULATION AREA AND BOUNDARY CONDITIONS . . . . . . . . . . . . . . . . . 5 4.3 NUMERICAL ANALYSYS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 4.4 FAST FOURIER TRANSFORMATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 4.5 INTERACTION OF FIELD AND PARTICLES . . . . . . . . . . . . . . . . . . . . . . . . 6 4.6 EVOLUTION OF FIELDS AND PARTICLES . . . . . . . . . . . . . . . . . . . . . . . . 7

4.6.1 predictors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 4.6.2 correctors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

5 INITIAL ELECTRON BEAM 9 6 INITIAL LASER FIELD 10 7 FREE SPACE FOR LASER FIELD 10 8 MIRROR REFRECTION 11 9 DRIFT SPACE FOR ELECTRON 12 10 DISPERSIVE SECTION 12 11 Appendix 13

11.1 [JJ] for m = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1 Abstract

A three-dimensional free electron laser simulation code, ELFIN, is introduced. The scheme to solve FEL equations and some result to show this

2 INTRODUCTION

ELFIN (Electron and Laser Field INteraction code) is a three-dimensional free electron laser simulation code for FEL interaction between a monochromatic laser and a continuous electron beam.

This is a paricle code [5].

The evolution of laser fields are described by a partial differential equation(PDE). In numeircal simulation, this PDE are transformed to a set of ordinally differential equations(ODE) by several methods.

One of the methods for it is mode-expansion ,in which laser fields are expanded with orthogonal modes such as TEM modes or laser cavity mode [1] and the PDE is transformed to a set of ODEs of the coefficients of the modes. Another method is a space discretization in the transverse space. A laser field is evaluated on the mesh points placed on the transvere space and PDE is transformed to a set of ODEs of the values of the laser field on the mesh points.

Mode-expansion is adequate for a cavity mode analysis but it is hard to apply to high-gain FELs in which strong optical guiding effects exist.

In space discretization scheme, the set of ODEs is usually solved by linear matrix solvers and an FEL with optical guiding effect can be treated.

But in this scheme, there exist the fast oscillation modes along the direction of the propagation of the laser field if the small spatial mesh size is used in transverse plane. Physically, this is a diffraction effect.

If we use large step size for the direction of the propagation with small mesh size in transverse plane, aliasing problems may cause unphysical field because of this fast oscillating modes.

The seed of this unphysical field mainly come from shot noise by worse statistics of superparticle than real electrons and it also unphysical.

To overcome this problem, Fourier mode expansion of a laser field with fast Fourier transformation is used in ELFIN to transform the PDE of a laser field to a set of ODEs of coefficients of Fourier modes, and these ODEs of the laser field and the ODEs for electrons are solved simultaneously with a kind of predictor-corrector methods. With this scheme, this code is almost free from aliasing problem as mentioned below.

For the propagation of laser field without electron beams, the result of the evolution of the laser field to the direction of its propagation is rigorous except that the transverse plane is discretized.

And it is easy to connect to the simulation of free space propagation of laser fields with the method shown in [4].

ELFIN consists of several modules shown in table 1. Each of them simulates a component of a free electron laser system such as FEL amplifiers, mirros, and free spaces.

SI unit is used in this paper. c is the speed of light and e, m_{e} are charge, mass and classical radius of electrons, respectively.

The coordinate system used here is as follows. The central axis of a wiggler is the z axis and the propagation of the laser field is the direction of z axis. The x axis is the direction of the wiggling of electrons. The y axis is set to form the left hand cartesian system in the relation with the x and z axis.

x̂, ŷ, ẑ are the unit vector in x, y, z direction, respectively.

To describe the state of electrons, the phase of the particles in a pondermotive potential made by the laser field of odd m-th harmonic, ζ = m ∫^{z} k_{W} dz + k_{L}z − ω_{L}t, the relativistic factor γ = E/m_{e}c^{2}, the transverse position and normalized transverse velocity, (x, y, β_{x} = v_{x}/c, β_{y} = v_{y}/c) are used.

3 ASSUMPTION

The following is the assumptions on which ELFIN is based.

table 1

Elenent  Function
FEL  Interaction of a lasers and electrons in a wiggler.
INITIAL LASER FIELD  Creation of an initial laser field.
FREE SPACE FOR LASER FIELD Propagation of a laser field in a free space.
MIRROR  Reflection of a laser field at a mirror.
ELECTRON BEAM  Creation of an electron beam.
DRIFT SPACE FOR ELECTRON Propagation of an electron beam in a drift space.
DISPERSIVE SECTION  Propagation of a laser field and an electron beam

in a dispersive section.

• An FEL system is time independent and z is used as the independent variable, which means that the laser field and the electron beam is static which means that the laser field is monochromatic and the electron beam is continuous.

But for an FEL with a bunched beam, if the bunch length of the electron beam and the pulse length of the laser field are much longer than N_{W} λ, where N_{W} is the period of a wiggler and λ is the wave length of the laser, we can treat it to be time-independent system locally and ELFIN can be adapted.

• The vector potential of a laser field has the form

A_{L}(x, y, z) = x̂ |A_{L}(x, y, z)| cos (k_{L}z − ω_{L}t + ϕ_{L}(x, y, z))  (1)
A_{L}(x, y, z) = |A_{L}(x, y, z)|e^{iϕ_{L}(x,y,z)}  (2)

.

The dimensionless amplitude of the vector potential of the laser field

a_{L}(x, y, z) = \frac{eA_{L}(x, y, z)}{\sqrt{2}m_{e}c}  (3)

is used to describe a status of the laser field.

• The vector potential of a wiggler field has the form of

A_{W} (x, y, z) = −x̂|A_{W} (x, y, z)| sin (∫ ^{z} k_{W} (z′)dz′ + ϕ_{W} (x, y, z))  (4)
A_{W} (x, y, z) = |A_{W} (x, y, z)|e^{iϕ_{W} (x,y,z)}.  (5)

. The dimensionless amplitude of vector potential of the wiggler is

a_{W} (x, y, z) = \frac{eA_{W} (x, y, z)}{\sqrt{2}m_{e}c}.  (6)

For static wiggler fields, ϕ_{W} can be set to 0 and the dependence of the field strength to the x, y co-ordinates is assumed to be

A_{W} (x, y, z) = A_{W 0}(z) (1 + \frac{1}{2}h_{x}(z)^{2}x^{2} + \frac{1}{2}h_{y}(z)^{2}y^{2})  (7)
h_{x}(z)^{2} + h_{y}(z)^{2} = k_{W} (z)^{2}.  (8)

• External magnetic fields has the form of

(B_{Ex}(x, y, z), B_{Ey}(x, y, z), B_{Ez}(x, y, z)) .  (9)

• External force has the form of

(F_{x}(x, y, z), F_{y}(x, y, z), F_{z}(x, y, z)) .  (10)

• The slowly varying amplitude and phase approximation is used. This approximaion is that the scale of the variying of the amplitude and the phase of the laser field are large enough compared with the wave length of the laser field;

\frac{λ_{L}}{|a_{L}|} \frac{∂|a_{L}|}{∂z}, \frac{λ_{L}}{|ϕ_{L}|} \frac{∂|ϕ_{L}|}{∂z} ≪ 1.  (11)

• ELFIN assume that the variables of a system does not change for the λ_{W} change of z, which means that for a quantitiy Q,

\frac{λ_{W}}{|Q|} \frac{∂|Q|}{∂z} ≪ 1.  (12)

.

With this assumption, ELFIN uses quantities averaged in the regard with wiggling motions of the electrons.

4 FEL AMPLIFIER

4.1 FEL EQUATIONS

The FEL equations used in ELFIN for m-th harmonics are

(−\frac{i}{2k_{L}} ∇_{⊥}^{2} + \frac{∂}{∂z} ) a_{L}(x, y, z) = s(P, x, y, z)  (13)

s(P, x, y, z) = ∑_{0≤ζ_{i}<2π} s_{P} (P_{i}, z)δ(x − x_{i})δ(y − y_{i})

s_{P} (P_{i}, z) = −r_{e} \frac{1}{γ_{i}} \{|a_{W i}|[JJ]_{i}e^{−i(ζ_{i}+ϕ_{Li})} + ia_{Li}\}  (14)

\frac{dγ_{i}}{dz} = k_{L} \frac{1}{γ_{i}} [JJ]_{i}|a_{W i}||a_{Li}| cos (ζ_{i} + ϕ_{W i} + ϕ_{Li}) + \frac{1}{m_{e}c^{2}} F_{z,i}  (15)

(16)

\frac{dζ_{i}}{dz} = mk_{W}  (17)

−k_{L} \frac{1}{2γ_{i}^{2}} \{1 + |a_{W i}|^{2} − 2|a_{W i}||a_{Li}|[JJ]_{i} sin (ζ_{i} + ϕ_{W i} + ϕ_{Li})\}

−\frac{k_{L}}{2} (β_{xi}^{2} + β_{yi}^{2})  (18)
\frac{dx_{i}}{dz} = β_{xi}  (19)
\frac{dy_{i}}{dz} = β_{yi}  (20)
\frac{d(γ_{i}β_{xi})}{dz} = −\frac{|a_{W 0}|^{2}h_{x}^{2}}{γ_{i}} x_{i} + \frac{e}{m_{e}c} (β_{yi}B_{Ez,i} − B_{Ey,i}) + \frac{1}{m_{e}c^{2}} F_{x,i}  (21)
\frac{d(γ_{i}β_{yi})}{dz} = −\frac{|a_{W 0}|^{2}h_{y}^{2}}{γ_{i}} y_{i} + \frac{e}{m_{e}c} (B_{Ex,i} − β_{xi}B_{Ez,i}) + \frac{1}{m_{e}c^{2}} F_{y,i}  (22)

(23) and

a_{Li} = a_{L}(x_{i}, y_{i}, z)  (24)
ϕ_{Li} = ϕ_{L}(x_{i}, y_{i}, z)  (25)
a_{W i} = a_{W} (x_{i}, y_{i}, z)  (26)
ϕ_{W i} = ϕ_{W} (x_{i}, y_{i}, z)  (27)
a_{W 0} = a_{W} (0, 0, z)  (28)
(B_{Ex,i}, B_{Ey,i}, B_{Ez,i}) = (B_{Ex}(x_{i}, y_{i}, z), B_{Ey}(x_{i}, y_{i}, z), B_{Ez}(x_{i}, y_{i}, z))  (29)
(F_{x,i}, F_{y,i}, F_{z,i}) = (F_{x}(x_{i}, y_{i}, z), F_{y}(x_{i}, y_{i}, z), F_{z}(x_{i}, y_{i}, z))  (30)
[JJ]_{i} = J_{\frac{m−1}{2}} (mξ_{i}) − J_{\frac{m+1}{2}} (mξ_{i})  (31)
ξ_{i} = \frac{|a_{W i}|^{2}}{2(1 + |a_{W i}|^{2})}  (32)
i = 1, 2, 3, …, N_{e}

, where J_{0}, J_{1} are the 0th order and 1st order Bessel functions, respectively. N_{e} is the number of electrons in the computation area and i = 1, 2, 3, …, N_{e}.

The equations (16)-(18) describe the longitudinal motions of the electrons and (19)-(23) describe the transverse motions of the electrons.

4.2 CALCULATION AREA AND BOUNDARY CONDITIONS

ELFIN analyzes the laser field and the source field in a square area of the size L_{a} ×L_{a} in transverse space (x, y). A periodic boundary condition

a_{L}(x, y, z) = a_{L}(x + L_{a}, y, z) = a_{L}(x, y + L_{a}, z)  (33)

is imposed on the laser field and the source field in this area.

Using this periodic boundary condition, we can check the validity of the size of the area.

If the value of the laser field and the value of the source field are both negligiblly small at the boundary compared with the central region, then the size of area is wide enough. If not, the area imposed is too small and the user should expand the area.

The laser field and the electron beam are static in time hence the periodic boundary condition is assumed for ζ coodinate of electrons with the period 2π.

4.3 NUMERICAL ANALYSYS

To solve these equations nummerically, ELFIN uses the Fourier transformation to transform the PDE (13) to a set of ODEs. Superparticles whose number is much smaller than the actual number of electrons is used.

ELFIN adopt the fast Fourier transformation(FFT) for this Fourier transformation. To apply FFT, the equidistant values of the fields are required hence the square mesh is imposed on (x, y) space and the laser field a_{L}(x, y, z) and the source field s(P, x, y, z) are evaluated on those mesh points.

The number of particle, N_{P}, are divided to N_{P} = N_{T} × N_{L} and the N_{T} position of the particles in the (x, y) space and N_{L} position of the particles in the γ space are randomly choosen with their assumed distribution by a Monte-Carlo method.

For ζ space, N_{L} positions of partcles are placed evenly in the region [0, 2π) to suppress unphysical spontaneous emission caused by artifitial randomness of the distribution.

The laser field and the particles interact each other through the Particle-In-Cell(PIC) method. These coupled ordinary differential equations of the field and the particles are solved with a kind of predictorcorrector mathods.

4.4 FAST FOURIER TRANSFORMATION

A equidistant square mesh is imposed on the area of size L_{a} × L_{a} in (x, y) space. The number of the mesh points is N_{a} × N_{a} and we set N_{a} is the power of two to apply FFT. The mesh size Δx is L_{a}/N_{a}.

The (x, y) coordinate of those mesh points are

(x_{p}, y_{q}) = (pΔx, qΔx)  (34)
p, q = −\frac{N_{a}}{2} + 1, −\frac{N_{a}}{2} + 2, …, −1, 0, 1, 2, …, \frac{N_{a}}{2}  (35)

In the following, p, q run in this range. The fields are evaluated on the mesh points (x_{p}, y_{q}). The Fourier modes in (x, y) space for the periodic boundary condition (33) are

u_{lm}(x, y) = exp (i\frac{2π}{L_{a}} (lx + my))  (36)

l, m = −\frac{N_{a}}{2} + 1, −\frac{N_{a}}{2} + 2, …, −1, 0, 1, 2, …, \frac{N_{a}}{2} .

In the following, l, m run in this range, l, m = −\frac{N_{a}}{2} + 1, −\frac{N_{a}}{2} + 2, …, −1, 0, 1, 2, …, \frac{N_{a}}{2} .

The Fourier coefficients of of the laser field, ã_{lm}(z) , and those of the source field, s̃_{lm}(z), are

ã_{lm}(z) = ∑_{p} ∑_{q} a_{L}(x_{p}, y_{q}, z)u_{lm}(x_{p}, y_{q})  (37)
s̃_{lm}(z) = ∑_{p} ∑_{q} s(x_{p}, y_{q}, z)u_{lm}(x_{p}, y_{q})  (38)

a_{L}(x_{p}, y_{q}, z) = \frac{1}{N_{a}^{2}} ∑_{l} ∑_{m} ã_{lm}(z)u_{lm}(x_{p}, y_{q})  (39)

s(x_{p}, y_{q}, z) = \frac{1}{N_{a}^{2}} ∑_{l} ∑_{m} s̃_{lm}(z)u_{lm}(x_{p}, y_{q}).  (40)

These transformation can be done with FFT by setting N_{a} = 2^{m}, where m is positive integer.

The partial differential equation of the laser field (13) is transformed by FFT to the set of ordinary differential equations of the Fourier coefficients ã_{lm}(z) and s̃_{lm}(z), as

(λ_{lm} + \frac{∂}{∂z} ) ã_{lm}(z) = s̃_{lm}(z)  (41)

λ_{lm} = \frac{i}{2k_{L}} (\frac{2π}{L_{a}} )^{2} (l^{2} + m^{2}).  (42)

The formal solution of them are

ã_{lm}(z + Δz) = e^{−λ_{lm}Δz} {ã_{lm}(z) + ∫_{z}^{z+Δz} s̃_{lm}(z′)e^{λ_{lm}z′}dz′} .  (43)

4.5 INTERACTION OF FIELD AND PARTICLES

N_{P} superparticles are used instead of N_{e} electrons and usually N_{e} >> N_{P}. The mass and charge of the particles are the N_{e}/N_{P} times larger than a electron and the equation of motion for these particles is the same as for electrons.

The Particle-In-Cell method(PIC) is used to interact the particles at the position of particles (x_{i}, y_{i}) and the laser field on the mesh points (x_{p}, y_{q}).

ELFIN interpolates the laser field a_{L}(x_{p}, y_{q}) on the mesh points to the positions of the particles (x_{i}, y_{i}) with the shape function W (x) as

a_{L}(x_{i}, y_{i}) = ∑_{p} ∑_{q} a_{L}(x_{p}, x_{q})W (x_{p} − x_{i})W (y_{q} − y_{i}).  (44)

The i-th particle creates the source field s_{i}(x, y, z) = s_{P} (P_{i}, z)δ(x − x_{i})δ(y − y_{i}). ELFIN replaces the δ functions in this source field with the shape function W (x), s_{i}(x, y, z) = s_{P} (P_{i}, z)W (x − x_{i})W (y − y_{i}) and the source field is evaluated on the mesh points (x_{p}, y_{q}) as

s(x_{p}, y_{q}, z) = ∑_{i=1}^{N_{P}} s_{P} (P_{i}, z)W (x_{p} − x_{i})W (y_{q} − y_{i}).  (45)

The second order spline function is used for the shape function and is

W (x) = \{ \begin{matrix}− (\frac{x}{Δx})^{2} + \frac{3}{4} & (|x| ≤ \frac{1}{2}Δx) \\ \frac{1}{2} (|\frac{x}{Δx}| − \frac{3}{2})^{2} & (\frac{1}{2}Δx < |x| ≤ \frac{3}{2}Δx) \\ 0 & (\frac{3}{2}Δx < |x|) \end{matrix}  (46)

4.6 EVOLUTION OF FIELDS AND PARTICLES

The position in z direction of each step is z_{n} = nΔz, n = 0, 1, 2, 3, …, N_{z}. where Δz is the step size of the z direction.

To evolve the laser field, we have to solve the equations (41), (λ_{lm} + \frac{∂}{∂z} ) ã_{lm}(z) = s̃_{lm}(z) ,which is the form of the equation of motion of forced oscillator and |\frac{1}{λ_{lm}} | is the oscillation period.

If we use naive methods to solve it, we have to use Δz shorter than δz = 1/|λ_{N_{a}N_{a}}| = 2k_{L}(\frac{Δx}{2π})^{2} for the transverse mesh size Δx. δz is the scale of the evolution of the laser field of the transverse scale Δx and is the order of the Rayleigh length.

However, The electron distribution in the transverse space has usually smooth function such as Gaussian or Lorentzian and its scale is the r.m.s. or HWHM of the distribution and are much larger than Δx. This means that the term s̃_{lm}(z) with large λ_{lm} have small values and is expected to be negligible.

The scale of change of the beam profile in z direction is the betatron period and the scale of the lontitudinal motion of the particles is synchrotron period. Hence ‘HLThese period have the scale of the length of the wiggler and large compared with δz. Based on these facts, Δz can be longer than δz if we can use as many particles as electrons.

But we use the particles whose number is much less than the real electrons and these particles have the transverse density fluctuation which is much larger than the real electrons. This fluctuation causes much larger s̃_{lm}(z) of large |λ_{lm}| than real situations. These s̃_{lm}(z) create ã_{lm}(z) of large |λ_{lm}|.

If one uses naive method to evolve the field with Δz ⪰ δz, aliasing problems may occurr for ã_{lm}(z) with large |λ_{lm}|.

To avoid the unphysical aliasing problem, ELFIN uses following method.

4.6.1 predictors

We will approximate the function s̃_{lm}(z) with the third-order polynomial which include the points (z_{n−3}, s̃_{lm}(z_{n−3}), (z_{n−2}, s̃_{lm}(z_{n−2}), (z_{n−1}, s̃_{lm}(z_{n−1})), (z_{n}, s̃_{lm}(z_{n})). This is

s̃_{lm}^{(p)}(z_{n} + uΔz) = ∑_{j=0}^{3} C_{j}^{(p)}u^{j}  (47)

where

C_{0}^{(p)} = s̃_{lm}(z_{n})  (48)
C_{1}^{(p)} = \frac{1}{6} (11s̃_{lm}(z_{n}) − 18s̃_{lm}(z_{n−1}) + 9s̃_{lm}(z_{n−2}) − 2s̃_{lm}(z_{n−3}))  (49)
C_{2}^{(p)} = \frac{1}{2} (2s̃_{lm}(z_{n}) − 5s̃_{lm}(z_{n−1}) + 4s̃_{lm}(z_{n−2}) − s̃_{lm}(z_{n−3}))  (50)
C_{3}^{(p)} = \frac{1}{6} (s̃_{lm}(z_{n}) − 3s̃_{lm}(z_{n−1}) + 3s̃_{lm}(z_{n−2}) − s̃_{lm}(z_{n−3})) .  (51)

We put this function into the integral in the equation (43) and integrate analytically from z_{n} to z_{n+1} = z_{n} + Δz, or u = 0 to 1. Then we get the approximated value of ã_{lm}(z) at z = z_{n+1}.

We call it the predictor of the laser field and express with

ã_{lm}^{(p)}(z_{n+1}) = exp (−λ_{lm}Δz) {ã_{lm}(z_{n}) + (∑_{j=0}^{3} C_{j}^{(p)}I_{lmj}) Δz}  (52)

where

I_{lmj} = ∫_{0}^{1} u^{j} exp (λ_{lm}uΔz) du.  (53)

This integrals are easily obtained using the equation

I_{lm,j+1} = \frac{∂I_{lm,j}}{∂ (λ_{lm}Δz)}  (54)

and are

I_{lm0} = (−1 + v) /w  (55)
I_{lm1} = (1 − v + vw) /w^{2}  (56)
I_{lm2} = (−2 + 2v − 2vw + vw^{2}) /w^{3}  (57)
I_{lm3} = (6 − 6v + 6vw − 3vw^{2} + vw^{3}) /w^{4}  (58)

(59)

where

v = exp (λ_{lm}Δz)  (60)
w = λ_{lm}Δz.  (61)

The predictor of the laser field, a_{L}^{(p)}(x, y, z_{n+1}), is obtained from ã_{lm}^{(p)}(z_{n+1}) by inverse Fourier transformation with FFT.

The equations of motion for particles (16)-(23) can be represented symbolically with the vecotr equation

\frac{dP(z)}{dz} = F (z, P(z), a_{L}(x, y, z))  (62)

where P(z) is

P(z) = (P_{1}(z), P_{2}(z), P_{3}(z), …, P_{N_{P}} (z))  (63)
P_{i}(z) = (ζ_{i}(z), γ_{i}(z), x_{i}(z), β_{xi}(z), y_{i}(z), β_{yi}(z))  (64)

and F (z, P(z), a_{L}(x, y, z)) is the left-hand-side term of the equation of motion (16)-(23). With this equation, We can obtain the predictor of the particles.

P^{(p)}(z_{n+1}) = P(z_{n}) + F (z_{n}, P(z_{n}), a_{L}(x, y, z_{n})) Δz.  (65)

Now we have the predictors of the laser field and the particles.

4.6.2 correctors

The corrector for the particles, P^{(c)}(z_{n+1}), which is more prousible value than predictor, can be obtain by the equation

P^{(c)}(z_{n+1}) = P(z_{n}) + \frac{1}{2} [F (z_{n}, P(z_{n}), a_{L}(x, y, z_{n})) + F (z_{n+1}, P^{(p)}(z_{n+1}), a_{L}^{(p)}(x, y, z_{n+1}))] Δz() (66)

and P^{(c)}(z_{n+1}) is taken as P(z_{n+1}).

Using this corrector of the particles P^{(c)}(z_{n+1}) and the predictor of the field a_{L}^{(p)}(x, y, z_{n+1}), we can get the Fourier coefficients of the source field s̃_{lm}^{(c)}(z_{n+1}).

Then we will approximate s̃_{lm}(z) with the third order polynomial which include the points, (z_{n−2}, s̃_{lm}(z_{n−2})), (z_{n−1}, s̃_{lm}(z_{n−1})), (z_{n}, s̃_{lm}(z_{n})), (z_{n+1}, s̃_{lm}^{(c)}(z_{n+1})) and it is

s̃_{lm}^{(c)}(z_{n} + uΔz) = ∑_{j=0}^{3} C_{j}^{(c)}u^{j}  (67)

where

C_{0}^{(c)} = s̃_{lm}(z_{n})  (68)
C_{1}^{(c)} = \frac{1}{6} (3s̃_{lm}(z_{n}) − 6s̃_{lm}(z_{n−1}) + s̃_{lm}(z_{n−2}) + 2s̃_{lm}^{(c)}(z_{n+1}))  (69)
C_{2}^{(c)} = \frac{1}{2} (−2s̃_{lm}(z_{n}) + s̃_{lm}(z_{n−1}) + s̃_{lm}^{(c)}(z_{n+1}))  (70)
C_{3}^{(c)} = \frac{1}{6} (−3s̃_{lm}(z_{n}) + 3s̃_{lm}(z_{n−1}) − s̃_{lm}(z_{n−2}) + s̃_{lm}^{(c)}(z_{n+1})) .  (71)

These s̃_{lm}^{(c)}(z) are put into the integral of the equation (43) and are analytically integrated, then we get the corrector,ã_{lm}^{(c)}(z_{n+1}), which is more plausible value than the predictor ã_{lm}^{(p)}(z_{n+1}),

ã_{lm}^{(c)}(z_{n+1}) = exp (−λ_{lm}Δz) {ã_{lm}(z_{n}) + (∑_{j=0}^{3} C_{j}^{(c)}I_{lmj}) Δz} .  (72)

We may continue this procudure to get more plousible value replacing this correctors of the particle and field by the predictor in (66). We set the obtained values by this method to be the values at the next step, z = z_{n+1}.

The laser field a_{L}(x, y, z_{n+1}) can be obtaned with the inverse Fourier transformation from ã_{lm}(z_{n+1}).

ELFIN can get the values ã_{lm}(z_{n+1}), P(z_{n+1}) and s̃_{lm}(z_{n+1}) from s̃_{lm}(z_{n−3}), s̃_{lm}(z_{n−2}), s̃_{lm}(z_{n−1}), s̃_{lm}(z_{n}), ã_{lm}(z_{n}) and P(z_{n}).

At z = z_{0} which is the initial point of the calculation, we do not have the history of s̃_{lm}(z). In this case, we set s̃_{lm}(z_{−3}) = s̃_{lm}(z_{0}) , s̃_{lm}(z_{−2}) = s̃_{lm}(z_{0}) and s̃_{lm}(z_{−1}) = s̃_{lm}(z_{0}) , assuming that step size we choose is small enough to do so.

With this predictor-corrector method, we can average the contribution s̃_{lm}(z) with large |λ_{lm}| hence the possibility of the unphysical aliasing problem is expected to be small with this method.

5 INITIAL ELECTRON BEAM

ELFIN uses the distribution of the form

D(ζ, γ, x, β_{x}, y, β_{y})dζ dγ dx dβ_{x} dy dβ_{y} = g(ζ)dζ h(γ)dγ f(x, β_{x}, y, β_{y})dx dβ_{x} dy dβ_{y}  (73)
f(x, β_{x}, y, β_{y})dx dβ_{x} dy dβ_{y} = N_{n} exp (−B^{n}) dx dβ_{x} dy dβ_{y},  (74)

B = (\frac{x − \overline{x}}{Δx})^{2} + (\frac{β_{x} − \overline{β_{x}}}{Δβ_{x}} )^{2} + (\frac{y − \overline{y}}{Δy} )^{2} + (\frac{β_{y} − \overline{β_{y}}}{Δβ_{y}} )^{2} ,  (75)

and g(ζ) is constant or uniform.

The particles which represent the electrons are generated with the Monte-Carlo method for the coodinate (γ, x, β_{x}, y, β_{y}).

The ζ coodinate of the particles are generated with the so-called ”quiet start” scheme. ,in which the particles are placed uniformly with equal spacing in ζ space. The fluctuation of the particle distrubution is much larger than the real situation because much smaller number of particles than the number of electrons. This flucutuation generate the unphysical spontaneous emission.

6 INITIAL LASER FIELD

The initial laser field is assumed to be TEM_{00} Gaussian mode with weist w_{0} ,the distance l_{f} to the weist and the power P_{0}.

The dimensionless vector potential of the laser field is

a_{L}(x, y, z) = \frac{a_{L0}}{1 + i\frac{z−l_{f}}{z_{R}}} exp (−\frac{1}{1 + i\frac{z−l_{f}}{z_{R}}} \frac{x^{2} + y^{2}}{w_{0}^{2}} )  (76)

(77)

,where λ_{L}z_{R} = πw_{0}^{2},

a_{L0} = \sqrt{P \frac{2r_{e}}{mc^{2}} \frac{1}{cπ^{2}}} \frac{λ_{L}}{w_{0}}
= 4.823382413 × 10^{−5}\sqrt{P [W]}\frac{λ_{L}}{w_{0}} .  (78)

z_{R} is Rayleigh length. The initial laser field on the mesh point is a_{L}(x_{p}, y_{q}, 0).

7 FREE SPACE FOR LASER FIELD

ELFIN can simulate the propagaion of the laser field in the free space of the different sizes of the calculation areas at the entrance and at the exit by Gaussian beam transformation [4].

The laser field propagating along with z axis in the free space obeys the equations

(−\frac{i}{2k_{L}} ∇_{⊥}^{2} + \frac{d}{dz} ) a_{L}(x, y, z) = 0  (79)

where L is the length of the free space and L_{in} and L_{out} are the sizes of calculation areas at the entrance and at the exit of the free space, respectively. The equation(79) is the quation(13) with s = 0. With Gaussian transformation,

( \begin{matrix}x′ \\ y′ \\ z′ − z′_{0} \end{matrix} ) = ( \begin{matrix}α\frac{x}{z+ib} \\ α\frac{y}{z+ib} \\ α^{2} \frac{z−z_{0}}{(z_{0}+ib)(z+ib)} \end{matrix} )  (80)

g(x, y, z) = \frac{1}{z + ib} exp (−ik_{L} \frac{x^{2} + y^{2}}{2(z + ib)})  (81)
v(x′, y′, z′) = \frac{u(x, y, z)}{g(x, y, z)}  (82)

(83)

here , α is the complex constant and b and z_{0} are the real constatnts, the equation (79) is transformed to

(−\frac{i}{2k_{L}} ∇_{⊥}^{2}′ + \frac{∂}{∂z′} ) v(x′, y′, z′) = 0.  (84)

In ELFIN, b and z′_{0} are set to be 0.

The z coodinates at the entrance and exit are set to z_{1} and z_{2} = z_{1} + L, respectively, and we set α = z_{1}. Hence z′_{1} is 0. With this assumption, the transformation is now has the form

( \begin{matrix}x′ \\ y′ \\ z′ \end{matrix} ) = ( \begin{matrix}α\frac{x}{z} \\ α\frac{y}{z} \\ α\frac{z_{1}}{z} (z − z_{1}) \end{matrix} )  (85)

g(x, y, z) = \frac{1}{z} exp (−ik_{L} \frac{x^{2} + y^{2}}{2z} )  (86)
v(x′, y′, z′) = \frac{u(x, y, z)}{g(x, y, z)}.  (87)

With this transformation, the free space with the width L_{in} and L_{out} of the areas at the entrance z = z_{0} and the exitz = z_{1}, respectively in the (x, y, z) coordinate system can be transformed to the free space of length L′ = \frac{z_{1}L}{z_{1}+L} and the same width L_{in} of the areas at the entrance and the exit in the (x′, y′, z′) coordinate system

To solve the equations (84), ELFIN uses the FFT meshod. First, the laser field on the mesh points (x′_{p}, y′_{q}) is Fourier transformed to their Fourier coefficients. The equations for them are

(λ_{lm} + \frac{∂}{∂z′} ) ã_{lm}(z′) = 0  (88)

and thier solutions are

ã_{lm}(z′ = L) = exp (−λ_{lm}L) ã_{lm}(z′ = 0)  (89)

(90)

The laser field a_{L}(x, y, z′ = L) can be obtained from ã_{lm}(z′ = L) by inverse Fourier transformation (87).

8 MIRROR REFRECTION

The parameters of a mirror which reflects a laser field are the radius r, curvature ρ, the (x, y) coordinate of the center of the mirror, , the transverse component of the normal vector θ_{x}, θ_{y} which represent the tilt of the mirror and the reflectance R. The contribution of the curvature of the mirror to the z coordinate of the surface of the mirror is

z_{ρ}(x, y) = \frac{r(x, y)^{2}}{2ρ}  (91)
r(x, y) = \sqrt{(x − Δx)^{2} + (y − Δy)^{2}}.

The contribution of the tilt is

z_{tilt}(x, y) = −xθ_{x} − yθ_{y}.  (92)

These amount are assumed to be small enough the z coordinate of the surface is just the sum of them and is

z(x, y) = z_{ρ}(x, y) + z_{tilt}(x, y)  (93)
= \frac{(x − Δx)^{2} + (y − Δy)^{2}}{2ρ} − xθ_{x} − yθ_{y}.  (94)

The phase difference Δ(x, y) of the laser field at point (x, y) is

Δ(x, y) = k_{L}d(x, y)
= −2k_{L} (\frac{(x − Δx)^{2} + (y − Δy)^{2}}{2ρ} − xθ_{x} − yθ_{y})  (95)

,which is caused by the difference of length of the path of the field ,d(x, y) = −2z(x, y).

The reflectivity R(x, y) is set to be zero outside of the mirror. But the sharp edge of the mirror may causes the Gibbs phenomena hence to avoid it, the smoothing of the edge with second order splining is performed.

The complex reflectivity which include the phase change and the absorption of the laser field is defined with the equation

C(x, y) = R(x, y) exp (iΔ(x, y))
= R(x, y) exp (2ik_{L} (−\frac{(x − Δx)^{2} + (y − Δy)^{2}}{2ρ} + xθ_{x} + yθ_{y}))  (96)

The reflected laser field a_{L}^{REF} is obtained with the equation

a_{L}^{REF} (x_{p}, y_{q}) = C(x_{p}, y_{q})a_{L}^{IN}(x_{p}, y_{q})  (97)

where a_{L}^{IN} is the incident laser field.

9 DRIFT SPACE FOR ELECTRON

The motion of electrons in a drift space of length L is

x^{out} = x^{in} + β_{x}^{in}L  (98)
y^{out} = y^{in} + β_{y}^{in}L  (99)

(100)

where superscript in and out show that the values with them are entrance and exit of this section, respectively.

10 DISPERSIVE SECTION

To make an FEL with an optical klystron configuration, A dispersive section is required between two fel sections.

The first FEL section produces the energy modulation of electrons.

The dispersive section, where the pass length of electrons depends on their energy, enhances longitudinal bunching of an electron beam.

The strength of the dispersive section is specified by how much phase of the laser field passes electrons in the dispersive section and the phase difference Δϕ between the entrance and at the exit of the dispersive section is used for it .

An dispersive section is specified with two parameters, an actual length and an effective length.

The actual legnth,L is the physical length of dispersive section.

The effective length, L_{d} is a length required to produce the phase difference Δϕ without any electric/magnetic field and L_{d} and Δζ have a relation,

Δζ = k_{L} (L_{d} − \frac{c}{v} L_{d}) = \frac{k_{L}L_{d}}{2γ^{2}}  (101)

where γ is the gamma factor of the electron beam, k_{L} is the wave number of the laser field, respectively.

In some papers, N_{d}, which is how many period the laser field passes electron in dispersive section, is used and is defined as Δζ = −2πN_{d}, hence we have

N_{d} = \frac{k_{L}L_{d}}{4πγ^{2}} .  (102)

There are several scheme to get longer effective length with shorter actual length. If the transverse magnetic field B_{y} used in a dispersive section to get longer effective length, then L_{d} is

L_{d} = L [1 + \frac{e^{2}}{m^{2}c^{2}L} ∫_{0}^{L} (∫_{0}^{z′} B_{y}(z′′)dz′′) dz′] .  (103)

In simulation, The laser field propagate this section just a free space of the length L and the change of parameters of the electron beam are

ζ^{out} = ζ^{in} − \frac{k_{L}L_{d}}{2γ^{2}} − \frac{k_{L}L}{2} (β_{x}^{in}^{2} + β_{y}^{in}^{2})  (104)
x^{out} = x^{in} + β_{x}^{in}L  (105)
y^{out} = y^{in} + β_{y}^{in}L  (106)

(107)

where superscript in and out show that the values with them are entrance and exit of this section, respectively. .

11 Appendix

11.1 [JJ] for m = 1

[JJ] for 1st harmonics is calculated through least square fitted 3rd-order polynomials for speed-up. The relative difference is less than 1.31 × 10^{−5}.

References

[1] For example, C. M. Tang and P. A. Sprangle, IEEE J. QUantum Electron. QE-21(1985),970.

[2] C. M. Tang, P. Sprangle, A. Ting and B. Hafizi , J.Appl.PHys. 66(4),15(1989)

[3] B. D. McVey, ”THREE-DIMENSIONAL SIMULATIONS OF FREE ELECTRON LASER PHYSICS,” Nucl. Instr. and Meth. A250(1986)449.

[4] E. A. Sziklas and A. E. Siegman, Proc. IEEE 62,410(1974).

[5] For example, R. W. Hockney and J. W. Eastwood, ”Computer Simulation using Particles,” Adam Higer. C. K. Birdsall and A. B. Langdon, ”Plasma Physics via Computer Simulation,” Adam Higer.