\documentstyle[titlepage,a4]{article}

\topmargin -1.5cm

\oddsidemargin  0.5cm
\evensidemargin 0.5cm

\textheight 23cm
\textwidth  16cm

\input{fel.inc}
\input{etc.inc}
\title{ELFIN\\
Three Dimensional Free Electron Laser Simulation Code}
\input{author.inc}

\date{\today}

\begin{document}

\maketitle
\tableofcontents
\newpage

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


\section{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 \cite{Hockney}.

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 \cite{Tang-Sprangle-85}
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 \cite{Sziklas}.

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.

$\vxe, \vye, \vze$ 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,
 $\zeta = \int^z k_W dz + k_L z - \omega_L t $, 
the relativistic factor $\gamma = E/m_e c^2$, 
the transverse position and normalized transverse velocity,
$\left( x, y, \beta_x = v_x / c,\beta_y = v_y / c \right)$ 
are used .

\begin{table}
\begin{center}
   table 1  \\
  \begin{tabular}{|l|l|}
    \hline
    Elenent & Function \\
    \hline
    \hline
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. \\
    \hline
  \end{tabular}
\end{center}
\end{table}

\section{ASSUMPTION}

The following is the assumptions on which ELFIN is based. 

\begin{itemize}

\item
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 \lambda$, where $N_W$ is the 
period of a wiggler and $\lambda$ is the wave length of the laser, 
we can treat it to be time-independent system locally 
and ELFIN can be adapted. 
 
\item
The vector potential of a laser field has the form 
\begin{eqnarray}
\vAL(x,y,z) &=& \vxe \left| A_L(x,y,z) \right| 
       \cos \left( k_L z-\omega_L t+\phi_L(x,y,z) \right) \\
A_L(x,y,z)  &=& |A_L(x,y,z)| e^{ i \phi_L(x,y,z) } 
\end{eqnarray}
. 

The dimensionless amplitude of the vector potential of the laser field 
\begin{equation} 
a_L(x,y,z)  =  \frac{e A_L(x,y,z)}{\sqrt{2} m_e c}   
\end{equation}
is used to describe a status of the laser field. 

\item 
The vector potential of a wiggler field has the form of 
\begin{eqnarray}
\vAW(x,y,z) & = &
 - \vxe |A_W(x,y,z)| 
    \sin\left(\dint{}{z}k_W(z')dz' + \phi_W(x,y,z) \right) \\
  A_W(x,y,z) & = & |A_W(x,y,z)| e^{ i \phi_W(x,y,z) }.  
\end{eqnarray}.
The dimensionless amplitude of vector potential of the wiggler is
\begin{equation} 
  a_W(x,y,z)  =  \frac{e A_W(x,y,z)}{\sqrt{2} m_e c} .
\end{equation}

For static wiggler fields, $\phi_W$ can be set to $0$ 
and the dependence of the field strength to the  $x,y$ co-ordinates is 
assumed to be
\begin{eqnarray}
A_W(x,y,z) &=& 
  A_{W0}(z)\left(1 + \half h_x(z)^2 x^2 + \half h_y(z)^2 y^2 \right) \\
   h_x(z)^2 + h_y(z)^2 & = & k_W(z)^2.
\end{eqnarray}

\item
External magnetic fields has the form of 
\begin{equation}
 \left(  B_{Ex}(x,y,z),B_{Ey}(x,y,z),B_{Ez}(x,y,z) \right). 
\end{equation}

\item
External force has the form of 
\begin{equation}
 \left(  F_x(x,y,z),F_y(x,y,z),F_z(x,y,z) \right). 
\end{equation}

\item
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;
\begin{eqnarray}
\frac{\lambda_L}{|a_L|}    \pdif{|a_L|}{z},
\frac{\lambda_L}{|\phi_L|} \pdif{|\phi_L|}{z}   & \ll &   1 .
\end{eqnarray}

\item
ELFIN assume that the variables of a system does not change
 for the $\lW$ change of z,
which means that for a quantitiy Q,
\begin{eqnarray}
    \frac{\lambda_W}{|Q|}  \pdif{|Q|}{z} & \ll &   1 .
\end{eqnarray}
.

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

The averaged value ${\overline Q}(z)$ for a quantity Q is 
\begin{equation}
{\overline Q}(z) = 
     \int_{z-\frac{\lambda_z}{2}}^{z+\frac{\lambda_z}{2}} Q(z') dz' 
\end{equation}
 
\end{itemize}

\section{FEL AMPLIFIER}

\subsection{FEL EQUATIONS}

The FEL equations used in ELFIN are 
\begin{equation}
   \left( -\frac{i}{2k_L}\nablat^2 + \pdif{}{z} \right) a_L(x,y,z)
=s(\vP,x,y,z)
                        \label{eq:mxwl}  
\end{equation}
\begin{eqnarray}
 s(\vP,x,y,z) & = & \sumzetai s_P(P_i,z) \delta (x-x_i ) \delta (y-y_i ) 
                        \nonumber \\
 s_P(P_i,z) & = &  -\re
 \frac{1}{\gamma_i} 
  \left\{ 
    |a_{Wi}| [JJ]_i e^{-i\left(\zeta_i+\phi_{Li}\right)}+ia_{Li} 
  \right\}
\end{eqnarray}

\begin{eqnarray} 
    \odif{\gamma_i}{z} & = & 
       k_L\frac{1}{\gamma_i} [JJ]_i |a_{Wi}| |a_{Li}|
          \cos\left( \zeta_i+\phi_{Wi}+\phi_{Li} \right) 
       +\frac{1}{m_e c^2} F_{z,i}\\
              \label{eq:eeqs} \label{eq:syns}\\
    \odif{\zeta_i}{z} & = & k_W \\
                      & & -  k_L \frac{1}{2\gamma_i^2} 
    \left\{ 
       1 + |a_{Wi}|^2 - 2|a_{Wi}||a_{Li}| [JJ]_i 
       \sin\left(\zeta_i + \phi_{Wi} + \phi_{Li} \right)
    \right\}   \nonumber \\
  &   & 
  - \frac{k_L}{2} \left( \beta_{xi}^2 + \beta_{yi}^2 \right) \label{eq:syne}\\ 
   \odif{x_i}{z} & = & \beta_{xi}  \label{eq:betas}\\
   \odif{y_i}{z} & = & \beta_{yi} \\
   \odif{(\gamma_i \beta_{xi})}{z} & = & 
             - \frac{|a_{W0}|^2 h_x^2}{\gamma_i} x_i 
             +\frac{e}{m_e c} \left( \beta_{yi} B_{Ez,i}-B_{Ey,i} \right)
             +\frac{1}{m_e c^2} F_{x,i} \\
   \odif{(\gamma_i \beta_{yi})}{z} &=& 
            - \frac{|a_{W0}|^2 h_y^2}{\gamma_i} y_i
             +\frac{e}{m_e c} \left( B_{Ex,i} - \beta_{xi} B_{Ez,i} \right) 
             +\frac{1}{m_e c^2} F_{y,i} \\
                   \label{eq:eeqe}  \label{eq:betae}
\end{eqnarray}
and
\begin{eqnarray}
   a_{Li}   & = & a_L (x_i,y_i,z)                 \\
   \phi_{Li} & = & \phi_L(x_i,y_i,z)                \\
   a_{Wi}   & = & a_W (x_i,y_i,z)                 \\
   \phi_{Wi} & = & \phi_W(x_i,y_i,z)                \\
   a_{W0}  & = & a_W(0,0,z)                        \\ 
   \left( B_{Ex,i},B_{Ey,i},B_{Ez,i} \right) &=&  
   \left( 
       B_{Ex}(x_i,y_i,z),B_{Ey}(x_i,y_i,z),B_{Ez}(x_i,y_i,z) 
   \right)\\ 
   \left( F_{x,i},F_{y,i},F_{z,i} \right) &=&  
   \left( 
       F_x(x_i,y_i,z),F_y(x_i,y_i,z),F_z(x_i,y_i,z) 
   \right)\\ 
   {[JJ]}_i  & = & J_0(\xi_i) - J_1(\xi_i)   \\
   \xi_i & = & \frac{|a_{Wi}|^2}{2(1+|a_{Wi}|^2)}    \\ 
       i  & = & 1,2,3,\ldots,N_e \nonumber 
\end{eqnarray}
, 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,\ldots,N_e$.

The equations \refeq{syns}-\refeq{syne} describe 
the longitudinal motions of the electrons 
and \refeq{betas}-\refeq{betae} describe the transverse motions
of the electrons.

\subsection{CALCULATION AREA AND BOUNDARY CONDITIONS}

ELFIN analyzes the laser field and the source field in a square area of 
the size $L_a \times L_a$ in transverse space $(x,y)$.
A periodic boundary condition 
\begin{equation}
   a_L(x,y,z) = a_L(x+L_a,y,z) =  a_L(x,y+L_a,z) 
      \label{eq:CyclicBoundaryCondition}
\end{equation}
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 $\zeta$ coodinate of 
electrons with the period $2\pi$.

\subsection{NUMERICAL ANALYSYS}

To solve these equations nummerically, 
ELFIN uses the Fourier transformation to transform 
the PDE \refeq{mxwl} 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(\vP,x,y,z)$ 
are evaluated on those mesh points.

The number of particle, $N_P$, are divided to $N_P = N_T \times N_L$ and
the $N_T$ position of the particles in the $(x,y)$ space
and $N_L$ position of the particles in the $\gamma$ space
are randomly choosen 
with their assumed distribution by a Monte-Carlo method.

For $\zeta$ space, 
$N_L$ positions of partcles 
are placed evenly in the region $[0,2\pi)$
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 predictor-corrector mathods. 


\subsection{FAST FOURIER TRANSFORMATION}

A equidistant square mesh is imposed on the area of 
size $L_a \times L_a$ in $(x,y)$ space.
The number of the mesh points is $N_a \times N_a $ 
and we set $N_a$ is the power of two to apply FFT.
The mesh size $\Delta x$ is $L_a/N_a$.

The $(x,y)$ coordinate of those mesh points are
\begin{equation}
 (x_p,y_q)  = \left( p \Delta x, q \Delta x \right)
\end{equation}
\begin{equation}
  p,q = 
  -\frac{N_a}{2}+1, -\frac{N_a}{2}+2,\ldots,-1,0,1,2,\ldots,\frac{N_a}{2}
\end{equation} 
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 \refeq{CyclicBoundaryCondition} are 
\begin{eqnarray}
 \ulm(x,y)  =  \exp \left( i \frac{2\pi}{L_a}(lx+my) \right) \\
 l,m = -\frac{N_a}{2}+1, -\frac{N_a}{2}+2,\ldots,-1,0,1,2,\ldots,\frac{N_a}{2} 
\nonumber.
\end{eqnarray}
In the following, $l,m$ run in this range,
$l,m = -\frac{N_a}{2}+1, -\frac{N_a}{2}+2,\ldots,-1,0,1,2,\ldots,\frac{N_a}{2}$. 


The Fourier coefficients of of the laser field, $\alm(z)$ , and those of  
 the source field, $\slm(z)$, are
\begin{eqnarray}
\alm(z) & = &   \sumpq a_L (x_p,y_q,z) \ulm(x_p,y_q) \\
\slm(z) & = &   \sumpq s(x_p,y_q,z)   \ulm(x_p,y_q) \\
a_L(x_p,y_q,z) & = & \frac{1}{N_a^2} \sumlm \alm(z) \ulm(x_p,y_q) \\
s(x_p,y_q,z)   & = & \frac{1}{N_a^2} \sumlm \slm(z) \ulm(x_p,y_q).
\end{eqnarray}
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 \refeq{mxwl} is 
transformed by FFT to the set of ordinary differential equations 
of the Fourier coefficients $\alm(z)$ and $\slm(z)$, 
as
\begin{equation}
   \left( \lambda_{lm} + \pdif{}{z} \right) \alm(z) = \slm(z) 
            \label{eq:almandlambda}
\end{equation}
\begin{equation}
   \lambda_{lm} = \frac{i}{2 k_L} \left(\frac{2\pi}{L_a}\right)^2 (l^2 +m^2 ).
       \label{eq:lambda_lm} 
\end{equation}
The formal solution of them are
\begin{equation}
\alm(z+\Delta z) =  
   e^{-\lambda_{lm} \Delta z } 
    \left\{
      \alm(z) 
        + \int_z^{z+\Delta z} \slm(z') e^{\lambda_{lm} z'} dz' 
    \right\}. 
\label{eq:sint}
\end{equation}

\subsection{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
\begin{equation}
a_L(x_i,y_i) = \sumpq a_L(x_p,x_q) W(x_p-x_i)W(y_q-y_i).
\end{equation}

The i-th particle creates the source field 
 $s_i(x,y,z) = s_P(P_i,z) \delta(x-x_i) \delta(y-y_i)$.
ELFIN replaces the $\delta$ 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 
\begin{eqnarray}
      s(x_p,y_q,z) &=& \sumip s_P(P_i,z) W(x_p-x_i) W(y_q-y_i). 
\end{eqnarray} 

The second order spline function is used for the shape function 
and is
\begin{equation}
W(x) = 
\left\{
   \begin{array}{ll}
         -\left( \frac{x}{\Delta x}\right)^{2} + \frac{3}{4} 
                                         & 
                                           ( |x| \leq \half\Delta x) \\
        \half \left( |\frac{x}{\Delta x}| - \frac{3}{2} \right)^2 
                                         & 
              \left( \half\Delta x < |x| \leq \frac{3}{2}\Delta x \right)\\
        0                                & 
              \left( \frac{3}{2}\Delta x < |x| \right) 
  \end{array}
\right. 
\end{equation}

\subsection{EVOLUTION OF FIELDS AND PARTICLES} 
 
The position in z direction of each step is 
$z_n = n \Delta z, n=0,1,2,3,\ldots, N_z$. 
where $\Delta z$ is the step size of the z direction.

To evolve the laser field,
we have to solve the equations \refeq{almandlambda},
$ \left( \lambda_{lm} + \pdif{}{z} \right) \alm(z) = \slm(z)$
,which is the form of the equation of motion of forced oscillator and
$|\frac{1}{\lambda_{lm}}|$ is the oscillation period.

If we use naive methods to solve it,
we have to use $\Delta z$ shorter than
$\delta z = 1/| \lambda_{N_a N_a} | = 2 k_L (\frac{\Delta x}{2\pi})^2 $
for the transverse mesh size $\Delta x$.
$\delta z$ is the scale of the evolution of the laser field 
of the transverse scale $\Delta 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 $\Delta x$.
This means that the term $\slm(z)$ with
large $\lambda_{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 $\delta z$.
Based on these facts, $\Delta z$ can be longer than $\delta 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 $\slm(z)$ of large $|\lambda_{lm}|$
 than real situations.
These $\slm(z)$ create $\alm(z)$ of large $|\lambda_{lm}|$.

If one uses naive method to evolve the field with $\Delta z \succeq \delta z$, 
aliasing problems may occurr for $\alm(z)$ with large $|\lambda_{lm}|$. 

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

\subsubsection{predictors}

We will approximate the function $\slm(z)$ with
the third-order polynomial which include the points 
$(z_{n-3},\slm({z_{n-3}})$, $(z_{n-2},\slm(z_{n-2})$, 
$(z_{n-1},\slm(z_{n-1}))$, $(z_n,\slm(z_n))$. 
This is
\begin{equation}
\slm^{(p)}(z_n + u\Delta z) = \sum_{j=0}^{3} C^{(p)}_{j}u^{j} 
\end{equation}
where
\begin{eqnarray}
C^{(p)}_{0}&=&\slm(z_n) \\
C^{(p)}_{1}
 &=& \frac{1}{6}
     \left(
          11\slm(z_n)-18\slm(z_{n-1})+9\slm(z_{n-2})-2\slm(z_{n-3}) 
     \right) \\
C^{(p)}_{2}
 &=& \frac{1}{2}
     \left(
           2\slm(z_n)- 5\slm(z_{n-1})+4\slm(z_{n-2})- \slm(z_{n-3})
     \right) \\
C^{(p)}_{3}
 &=& \frac{1}{6}
     \left(
            \slm(z_n)- 3\slm(z_{n-1})+3\slm(z_{n-2})- \slm(z_{n-3})
     \right).
\end{eqnarray} 
We put this function into the integral in the equation \refeq{sint}
and integrate analytically from $z_n$ to $z_{n+1}=z_n + \Delta z$, 
or $u = 0$ to $1$. 
Then we get the approximated value of $\alm(z)$ at $z = z_{n+1}$.

We call it the predictor of the laser field 
and express with 
\begin{equation}
  \alm^{(p)}(z_{n+1}) = \exp\left( -\lambda_{lm} \Delta z \right)
  \left\{ 
      \alm(z_n)+ \left( \sum_{j=0}^3 C^{(p)}_j I_{lmj} \right) \Delta z
  \right\}
\end{equation}
where
\begin{eqnarray}
 I_{lmj} & = & \dint{0}{1} u^j \exp\left(\lambda_{lm} u \Delta z\right) du .
\end{eqnarray}
This integrals are easily obtained using 
  the equation
\begin{equation}
  I_{lm,j+1} = \pdif{ I_{lm,j}}{\left(\lambda_{lm} \Delta z\right)}
\end{equation}
and are
\begin{eqnarray}
 I_{lm0} & = & \left( -1 +  v \right)/w                                    \\
 I_{lm1} & = & \left(  1 -  v + vw \right)/w^2                             \\
 I_{lm2} & = & \left( -2 + 2v - 2vw + vw^2 \right)/w^3                     \\
 I_{lm3} & = & \left(  6 - 6v + 6vw - 3vw^2 + vw^3 \right)/w^4           \\
\end{eqnarray}
where
\begin{eqnarray}
             v & = & \exp\left(\lambda_{lm} \Delta z \right)           \\
             w & = & \lambda_{lm} \Delta z .
\end{eqnarray}

The predictor of the laser field, $a_L^{(p)}(x,y,z_{n+1})$,
is obtained from $\alm^{(p)}(z_{n+1})$ by 
 inverse Fourier transformation with FFT.

The equations of motion for particles \refeq{eeqs}-\refeq{eeqe}
can be represented symbolically with the vecotr equation
\begin{equation}
   \frac{d\vP(z)}{dz}  =  \vF\left(z,\vP(z),a_L(x,y,z)\right)
\end{equation}
where $\vP(z)$ is  
\begin{eqnarray}
   \vP(z) = \left(
               P_1(z), P_2(z), P_3(z), \ldots, P_{N_P}(z) 
            \right) \\
   P_i(z) = \left(
               \zeta_i(z), \gamma_i(z)
               , x_i(z), \beta_{xi}(z), y_i(z), \beta_{yi}(z) 
            \right)
\end{eqnarray}
and $\vF\left(z,\vP(z),a_L(x,y,z)\right)$ is 
the left-hand-side term of the equation of motion \refeq{eeqs}-\refeq{eeqe}.
With this equation,
We can obtain the predictor of the particles.
\begin{equation}
 \vP^{(p)}(z_{n+1}) = \vP(z_n) 
         + \vF\left(z_n,\vP(z_n),a_L(x,y,z_n)\right) \Delta z .
\end{equation}

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

\subsubsection{correctors}

The corrector for the particles, $\vP^{(c)}(z_{n+1})$,
which is more prousible value than predictor, 
can be obtain by the equation 
\begin{equation}
\vP^{(c)}(z_{n+1}) = \vP(z_n) 
+ \frac{1}{2}
  \left[
       \vF\left(z_n,\vP(z_n),a_L(x,y,z_n)\right) 
          +\vF\left(z_{n+1},\vP^{(p)}(z_{n+1}),a_L^{(p)}(x,y,z_{n+1})\right) 
   \right] \Delta z
  \eq{corrector} 
\end{equation}     
and $\vP^{(c)}(z_{n+1})$ is taken as $\vP(z_{n+1})$.

Using this corrector of the particles $\vP^{(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 $\slm^{(c)}(z_{n+1})$.

Then we will approximate $\slm(z)$ with the third order polynomial
which include the points, 
 $(z_{n-2},\slm(z_{n-2}))$,
 $(z_{n-1},\slm(z_{n-1}))$,
 $(z_n,\slm(z_n))$, 
 $(z_{n+1},\slm^{(c)}(z_{n+1}))$
and it is 
\begin{equation}
\slm^{(c)}(z_n + u \Delta z) = \sum_{j=0}^{3} C^{(c)}_j u^j 
\end{equation}
where
\begin{eqnarray}
C^{(c)}_0 &=&   \slm(z_n) \\
C^{(c)}_1 &=& \frac{1}{6}
      \left(
         3 \slm(z_n) - 6\slm(z_{n-1}) + \slm(z_{n-2}) + 2\slm^{(c)}(z_{n+1})
      \right)   \\
C^{(c)}_2 &=& \frac{1}{2}
      \left(
        -2 \slm(z_n) + \slm(z_{n-1})                   +  \slm^{(c)}(z_{n+1})
      \right)  \\
C^{(c)}_3 &=& \frac{1}{6}
      \left(
        -3 \slm(z_n) + 3\slm(z_{n-1}) - \slm(z_{n-2})  +  \slm^{(c)}(z_{n+1})
      \right).
\end{eqnarray}

These $\slm^{(c)}(z)$ are put into the integral of the equation \refeq{sint}
and are analytically integrated, 
then we get the corrector,$\alm^{(c)}(z_{n+1})$,
which is more plausible value than the predictor $\alm^{(p)}(z_{n+1})$,
\begin{equation}
 \alm^{(c)}(z_{n+1}) = 
\exp\left( -\lambda_{lm} \Delta z \right)
 \left\{ 
   \alm(z_n) +  \left( \sum_{j=0}^3 C^{(c)}_j I_{lmj} \right) \Delta z
 \right\} .
\end{equation}

We may continue this procudure to get more plousible value 
 replacing this correctors of the particle and field by the predictor in \refeq{corrector}.
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 $\alm(z_{n+1})$.
  
ELFIN can get the values 
$\alm(z_{n+1})$, $\vP(z_{n+1})$ and $\slm(z_{n+1})$
from 
$\slm(z_{n-3})$, $\slm(z_{n-2})$, $\slm(z_{n-1})$, $\slm(z_n)$, $\alm(z_n)$
and $\vP(z_n)$.

At $z=z_0$ which is the initial point of the calculation,
we do not have the history of $\slm(z)$.
In this case,
 we set 
 $\slm(z_{-3}) = \slm(z_0)~$ , $\slm(z_{-2}) = \slm(z_0)~$ and $\slm(z_{-1}) = \slm(z_0)~$ ,
assuming that 
step size we choose is small enough to do so. 

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

\section{INITIAL ELECTRON BEAM}

ELFIN uses the distribution of the form 
\begin{equation}
  D(\zeta,\gamma,x,\beta_x,y,\beta_y)
      d\zeta \, d\gamma \,  dx \, d\beta_x \, dy \, d\beta_y
 = g(\zeta) d\zeta \; h(\gamma) d\gamma \; 
        f(x,\beta_x,y,\beta_y)dx \, d\beta_x \, dy \, d\beta_y
\end{equation}
\begin{equation}
 f(x,\beta_x,y,\beta_y) dx \, d\beta_x \, dy \, d\beta_y
       = N_n \exp\left(-B^n\right) dx \, d\beta_x \, dy \, d\beta_y,
\end{equation}

\begin{equation}
 B =  \left( \frac{x-\overline{x}}              {\Delta x} \right)^2 
    + \left( \frac{\beta_x-\overline{\beta_x} } {\Delta \beta_x} \right)^2 
    + \left( \frac{y-\overline{y}}              {\Delta y} \right)^2 
    + \left( \frac{\beta_y-\overline{\beta_y} } {\Delta \beta_y} \right)^2
         \nonumber ,
\end{equation}
and $g(\zeta)$ is constant or uniform.

The particles which represent the electrons are generated 
with the Monte-Carlo method for the coodinate 
$(\gamma,x,\beta_x,y,\beta_y)$. 

The $\zeta$ coodinate of the particles are generated 
with the so-called "quiet start" scheme. 
,in which
the particles are placed uniformly with equal spacing in $\zeta$ 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. 

\section{INITIAL LASER FIELD}

The initial laser field is assumed to be ${\rm 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 
\begin{eqnarray}
 a_L(x,y,z) & = &
   \frac{a_{L0}}{ 1+ i\frac{z-l_f}{z_R} }
       \exp\left( 
             -\frac{1}{ 1+ i\frac{z-l_f}{z_R} }\frac{x^2+y^2}{w_0^2} 
           \right)  \\
\end{eqnarray}
,where $\lambda_L z_R = \pi w_0^2$,
\begin{eqnarray}
  a_{L0} & = &
\sqrt{ P\frac{2\re}{mc^2 } \frac{1}{c\pi^2} } \frac{\lambda_L}{w_0}
                                                     \nonumber \\
 & = & 4.823382413 \times 10^{-5} \sqrt{P[W]} \frac{\lambda_L}{w_0}.
\end{eqnarray}
$z_R$ is Rayleigh length. 
The initial laser field on the mesh point is $a_L(x_p,y_q,0)$. 

\section{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 \cite{Sziklas}.

The laser field propagating along with z axis in the free space 
obeys the equations
\begin{equation}
  \left(
     -\frac{i}{2 k_L}\nablat^{2} + \frac{d}{dz} 
  \right) a_L(x,y,z) = 0 \label{eq:mxwl-free}
\end{equation}
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\refeq{mxwl-free} is the quation\refeq{mxwl} with $s=0$.
With Gaussian transformation,
\begin{eqnarray}
 \left(
    \begin{array}{c}
      x' \\
      y' \\
      z'-z_0'
    \end{array} 
 \right)
 &      =  &
 \left( 
    \begin{array}{c}
        \alpha \frac{x}{z+ib} \\
        \alpha \frac{y}{z+ib} \\
        \alpha^2 \frac{z-z_0}{(z_0+ib)(z+ib)}
    \end{array} 
 \right) \\
 g(x,y,z)  & = & \frac{1}{z+ib} 
        \exp\left( -i k_L \frac{x^2 + y^2 }{2 (z+ib)} \right)\\
 v(x',y',z') & = & \frac{u(x,y,z)}{g(x,y,z)}  \\
\end{eqnarray}
here , $\alpha$ is the complex constant
 and $b$ and  $z_0$ are the real constatnts,
the equation \refeq{mxwl-free} is transformed to
\begin{eqnarray}
\left( -\frac{i}{2 k_L}\nablat^2{}' + \pdif{}{z'} \right) v(x',y',z') & = & 0 
      \label{eq:mxwl-free-prime}.
\end{eqnarray}

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 $\alpha = z_1$.
Hence $z_1'$ is $0$.
With this assumption, the transformation is now has the form
\begin{eqnarray}
 \left( 
   \begin{array}{c}
      x' \\
      y' \\
      z'
   \end{array} 
 \right)
 & = &
 \left(
    \begin{array}{c}
       \alpha \frac{x}{z} \\
       \alpha \frac{y}{z} \\
       \alpha \frac{z_1}{z} ( z - z_1 )
    \end{array} 
 \right) \\
 g(x,y,z)  & = &  
    \frac{1}{z} \exp\left( -i k_L\frac{x^2 + y^2 }{2 z} \right) 
          \label{eq:g-sol} \\
 v(x',y',z')  & = & \frac{u(x,y,z)}{g(x,y,z)} \label{eq:g-trans}.
\end{eqnarray}
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 exit$z=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 \refeq{mxwl-free-prime},
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 
\begin{equation}
   \left( \lambda_{lm} + \pdif{}{z'} \right) \alm(z') = 0
\end{equation}
and 
thier solutions are
\begin{eqnarray}
   \alm(z'=L) & = & \exp\left( -\lambda_{lm} L \right) \alm(z'=0) \\
\end{eqnarray}

The laser field $a_L(x,y,z'=L)$ can be obtained from $\alm(z'=L)$ 
by inverse Fourier transformation \refeq{g-trans}. 

\section{MIRROR REFRECTION}

The parameters of a mirror which reflects a laser field 
are the radius $r$, curvature $\rho$, 
the $(x,y)$ coordinate of the center of the mirror,
, the transverse component of the normal vector  $\theta_x, \theta_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 
\begin{eqnarray}
 z_{\rho}(x,y) & = & \frac{r(x,y)^{2}}{2\rho}  \\
 r(x,y) & = & \sqrt{ (x-\Delta x)^{2} + (y-\Delta y)^{2} } .\nonumber
\end{eqnarray}
The contribution of the tilt is 
\begin{equation}
  z_{tilt}(x,y) = -x \theta_x - y \theta_y.
\end{equation}
These amount are assumed to be small enough  
the z coordinate of the surface is just the sum of them
and is
\begin{eqnarray}
  z(x,y) & = & z_{\rho}(x,y) + z_{tilt}(x,y) \\
  & = &   \frac{(x-\Delta x)^2 + (y-\Delta y)^2 }{2\rho}  
            -x \theta_x -y \theta_y.
\end{eqnarray}

The phase difference $\Delta (x,y)$ of the laser field at point $(x,y)$
is
\begin{eqnarray}
  \Delta(x,y) & = & k_L d(x,y) \nonumber \\ 
 & = & -2 k_L 
     \left(
       \frac{(x-\Delta x)^2 + (y-\Delta y)^2 }{2\rho} -x \theta_x -y \theta_y 
     \right)
\end{eqnarray}
,which is caused by the difference of length of the path of the field
,$d(x,y) = -2 z(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
\begin{eqnarray}
 C(x,y) & = &  R(x,y) \exp\left(i\Delta (x,y)\right) \nonumber \\
        & = &                               
 R(x,y)
 \exp\left( 2i k_L
         \left(
             -\frac{(x-\Delta x)^{2}+(y-\Delta y)^{2}}{2\rho} 
                  + x\theta_x + y\theta_y 
         \right)
      \right) 
\end{eqnarray}

The reflected laser field $a_L^{REF}$ is obtained with the equation
\begin{eqnarray}
 a_L^{REF}(x_p,y_q) & = & C(x_p,y_q) a_L^{IN}(x_p,y_q)  
\end{eqnarray}
where $a_L^{IN}$ is the incident laser field.

\section{DRIFT SPACE FOR ELECTRON}
The motion of electrons in a drift space of length $L$ is
\begin{eqnarray} 
    x^{out} &=& x^{in} + \beta_x^{in} L \\ 
    y^{out} &=& y^{in} + \beta_y^{in} L \\ 
\end{eqnarray}
where superscript in and out show that the values with them are
entrance and exit of this section, respectively.
 
\section{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 $\Delta\phi$
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 $\Delta\phi$ 
without any electric/magnetic field
and 
$L_d$ and $\Delta\zeta$ have a relation, 
\begin{equation}
   \Delta\zeta = k_L \left( L_d - \frac{c}{v} L_d  \right) 
               = \frac{k_L L_d}{2\gamma^{2}} 
\end{equation}
where $\gamma$ 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 
$\Delta \zeta = - 2\pi  N_d$, 
hence we have
\begin{eqnarray}
   N_d  &=& \frac{k_L L_d}{4\pi\gamma^{2}}  .
\end{eqnarray}

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
\begin{eqnarray}
   L_d &=& L 
       \left[
         1+\frac{e^2}{m^2 c^2 L}
            \int_{0}^{L}
               \left(
                  \int_{0}^{z'} B_y(z'') dz'' 
               \right) dz'
       \right]  .
\end{eqnarray}

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
\begin{eqnarray}
    \zeta^{out} &=& \zeta^{in}
        - \frac{k_L L_d}{2 \gamma^2} 
        - \frac{k_L L}{2}\left( {\beta_x^{in}}^2 + {\beta_y^{in}}^2 \right) \\
    x^{out} &=& x^{in} + \beta_x^{in} L  \\
    y^{out} &=& y^{in} + \beta_y^{in} L  \\
\end{eqnarray}
where superscript in and out show that the values with them are
entrance and exit of this section, respectively.
.

\thebibliography{99}
\bibitem{Tang-Sprangle-85}
For example,
C. M. Tang and P. A. Sprangle,
IEEE J. QUantum Electron. QE-21(1985),970.
\bibitem{Tang-Sprangle-89}
C. M. Tang, P. Sprangle, A. Ting and B. Hafizi ,
J.Appl.PHys. 66(4),15(1989)
\bibitem{McVey}
B. D. McVey,
"THREE-DIMENSIONAL SIMULATIONS OF FREE ELECTRON LASER PHYSICS,"
Nucl. Instr. and Meth. A250(1986)449.
\bibitem{Sziklas}
E. A. Sziklas and A. E. Siegman,
Proc. IEEE 62,410(1974).
\bibitem{Hockney}
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.

\end{document}
