1.5cm 1.5cm -0.75in 12pt 25pt 12pt 75pt 25cm 15cm \*(DY .nr PO 0 .AM .pl 100i .nr LL 7i .nr HF 0i .nr FM 0i .ds CH .ds CF .na .TL ELFIN .br Three Dimensional Free Electron Laser Simulation Code .AU Takeshi Nakamura .br SPring-8 .br Kamigori-cho, Ako-gun, Hyogo 678-12 JAPAN .br nakamura@spring8.or.jp .NH 1 Abstract .XS Abstract .XE .PP A three-dimensional free electron laser simulation code, ELFIN, is introduced. The scheme to solve FEL equations and some result to show this .NH 1 INTRODUCTION .XS INTRODUCTION .XE .PP ELFIN (Electron and Laser Field INteraction code) is a three-dimensional free electron laser simulation code and simulates a three-dimensional free electron laser(FEL) interaction between a monochromatic laser and a continuous electron beam. There are many three-dimensional FEL simulation codes in the world. Paricle methodsare used in most codes and ELFIN also adopt this method. The evolution of laser fields are described by a partial differential equation(PDE). In numeircal simulation, this PDE are converted to a set of ordinally differential equations(ODE) by several methods. One method to do it is mode-expansion method. In this method, laser fields are expanded with orthogonal modes such as TEM modes or cavity mode and the PDE is converted to a set of ODEs of the coefficients of the modes. Another method is to use the spatial mesh in the transverse space. A laser field is evaluated on the mesh points and PDE is converted to a set of ODEs of the values of the laser field on the mesh points. Formar scheme is good for cavity mode analysis but it is hard to apply to high-gain FELs in which strong optical guiding effects exist. In later 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 of the laser field 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 convert 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 . 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. .B x .R , .B y .R , .B z .R 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, = ^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 . .ce 100 table 1 .br .TS tab(&); L L . _ 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. _ .TE .ce 0 .NH 1 ASSUMPTION .XS ASSUMPTION .XE .PP The following is the assumptions on which ELFIN is based. .RS .IP \(bu 2 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. .IP \(bu 2 The vector potential of a laser field has the form .B A .R _L(x,y,z)= .B x .R | A_L(x,y,z) | ( k_L z-_L t+_L(x,y,z) ) .br A_L(x,y,z)= |A_L(x,y,z)| e^ i _L(x,y,z) . The dimensionless amplitude of the vector potential of the laser field a_L(x,y,z) = e A_L(x,y,z)2 m_e c is used to describe a status of the laser field. .IP \(bu 2 The vector potential of a wiggler field has the form of .B A .R _W(x,y,z) = - .B x .R |A_W(x,y,z)| (_^zk_W(z')dz' + _W(x,y,z) ) .br A_W(x,y,z) = |A_W(x,y,z)| e^ i _W(x,y,z). . The dimensionless amplitude of vector potential of the wiggler is a_W(x,y,z) = e A_W(x,y,z)2 m_e c . 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_W0(z)(1 + 12h_x(z)^2 x^2 + 12h_y(z)^2 y^2 ) .br h_x(z)^2 + h_y(z)^2 = k_W(z)^2. .IP \(bu 2 External magnetic fields has the form of ( B_Ex(x,y,z),B_Ey(x,y,z),B_Ez(x,y,z) ). .IP \(bu 2 External force has the form of ( F_x(x,y,z),F_y(x,y,z),F_z(x,y,z) ). .IP \(bu 2 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; _L|a_L| |a_L|z, _L|_L| |_L|z 1 . .IP \(bu 2 ELFIN assume that the variables of a system does not change for the _W change of z, which means that for a quantitiy Q, _W|Q| |Q|z 1 . . With this assumption, ELFIN uses quantities averaged in the regard with wiggling motions of the electrons. The averaged value Q(z) for a quantity Q is Q(z) = _z-_z2^z+_z2 Q(z') dz' .RE .NH 1 FEL AMPLIFIER .XS FEL AMPLIFIER .XE .PP .NH 2 FEL EQUATIONS .PP The FEL equations used in ELFIN are ( -i2k_L_^2 + z ) a_L(x,y,z) =s( .B P .R ,x,y,z) s( .B P .R ,x,y,z) = _0 _i< 2s_P(P_i,z) (x-x_i ) (y-y_i ) .br s_P(P_i,z) = -r_e 1_i |a_Wi| [JJ]_i e^-i(_i+_Li)+ia_Li d _id z = k_L1_i [JJ]_i |a_Wi| |a_Li| ( _i+_Wi+_Li ) +1m_e c^2 F_z,i .br .br d _id z = k_W .br - k_L 12_i^2 1 + |a_Wi|^2 - 2|a_Wi||a_Li| [JJ]_i (_i + _Wi + _Li ) .br + 12( _xi^2 + _yi^2 ) .br d x_id z = _xi .br d y_id z = _yi .br d (_i _xi)d z = - |a_W0|^2 h_x^2_i x_i +em_e c ( _yi B_Ez,i-B_Ey,i ) +1m_e c^2 F_x,i .br d (_i _yi)d z= - |a_W0|^2 h_y^2_i y_i +em_e c ( B_Ex,i - _xi B_Ez,i ) +1m_e c^2 F_y,i .br and a_Li = a_L (x_i,y_i,z) .br _Li = _L(x_i,y_i,z) .br a_Wi = a_W (x_i,y_i,z) .br _Wi = _W(x_i,y_i,z) .br a_W0 = a_W(0,0,z) .br ( 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) ) .br ( 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) ) .br [JJ]_i = J_0(_i) - J_1(_i) .br _i = |a_Wi|^22(1+|a_Wi|^2) .br 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 ()-() describe the longitudinal motions of the electrons and ()-() describe the transverse motions of the electrons. .NH 2 CALCULATION AREA AND BOUNDARY CONDITIONS .PP 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) 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. .NH 2 NUMERICAL ANALYSYS .PP To solve these equations nummerically, ELFIN uses the Fourier transformation to transform the PDE () 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( .B P .R ,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 predictor-corrector mathods. .NH 2 FAST FOURIER TRANSFORMATION .PP 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 ) p,q = -N_a2+1, -N_a2+2,...,-1,0,1,2,...,N_a2 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 () are u_lm(x,y) = ( i 2L_a(lx+my) ) .br l,m = -N_a2+1, -N_a2+2,...,-1,0,1,2,...,N_a2 . In the following, l,m run in this range, l,m = -N_a2+1, -N_a2+2,...,-1,0,1,2,...,N_a2. The Fourier coefficients of of the laser field, a_lm(z) , and those of the source field, s_lm(z), are a_lm(z) = _p _qa_L (x_p,y_q,z) u_lm(x_p,y_q) .br s_lm(z) = _p _qs(x_p,y_q,z) u_lm(x_p,y_q) .br a_L(x_p,y_q,z) = 1N_a^2 _l _ma_lm(z) u_lm(x_p,y_q) .br s(x_p,y_q,z) = 1N_a^2 _l _ms_lm(z) u_lm(x_p,y_q). 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 () is transformed by FFT to the set of ordinary differential equations of the Fourier coefficients a_lm(z) and s_lm(z), as ( _lm + z ) a_lm(z) = s_lm(z) _lm = i2 k_L (2L_a)^2 (l^2 +m^2 ). The formal solution of them are a_lm(z+z) = e^-_lm z a_lm(z) + _z^z+z s_lm(z') e^_lm z' dz' . .NH 2 INTERACTION OF FIELD AND PARTICLES .PP 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 _qa_L(x_p,x_q) W(x_p-x_i)W(y_q-y_i). 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_Ps_P(P_i,z) W(x_p-x_i) W(y_q-y_i). The second order spline function is used for the shape function and is W(x) = ll -( xx)^2 + 34 ( |x| 12x) .br 12( |xx| - 32 )^2 ( 12x < |x| 32x ) .br 0 ( 32x < |x| ) . .NH 2 EVOLUTION OF FIELDS AND PARTICLES .PP 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 (), ( _lm + z ) a_lm(z) = s_lm(z) ,which is the form of the equation of motion of forced oscillator and |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 | = 2 k_L (x2)^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 a_lm(z) of large |_lm|. If one uses naive method to evolve the field with z z, aliasing problems may occurr for a_lm(z) with large |_lm|. To avoid the unphysical aliasing problem, ELFIN uses following method. 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 + uz) = _j=0^3 C^(p)_ju^j C^(p)_0=s_lm(z_n) .br C^(p)_1 = 16 ( 11s_lm(z_n)-18s_lm(z_n-1)+9s_lm(z_n-2)-2s_lm(z_n-3) ) .br C^(p)_2 = 12 ( 2s_lm(z_n)- 5s_lm(z_n-1)+4s_lm(z_n-2)- s_lm(z_n-3) ) .br C^(p)_3 = 16 ( s_lm(z_n)- 3s_lm(z_n-1)+3s_lm(z_n-2)- s_lm(z_n-3) ) . We put this function into the integral in the equation () 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 a_lm(z) at z = z_n+1. We call it the predictor of the laser field and express with a_lm^(p)(z_n+1) = ( -_lm z ) a_lm(z_n)+ ( _j=0^3 C^(p)_j I_lmj ) z I_lmj = _0^1 u^j (_lm u z) du .br I_lm0 = ( -1 + v )/w .br I_lm1 = ( 1 - v + vw )/w^2 .br I_lm2 = ( -2 + 2v - 2vw + vw^2 )/w^3 .br I_lm3 = ( 6 - 6v + 6vw - 3vw^2 + vw^3 )/w^4 .br v = (_lm z ) .br w = _lm z . The predictor of the laser field, a_L^(p)(x,y,z_n+1), can be btained from a_lm^(p)(z_n+1) by inverse Fourier transformation with FFT. The equations of motion for particles ()-() can be represented with the equation d .B P .R (z) .R dz = .B F .R (z, .B P .R (z),a_L(x,y,z)) where .B P .R (z) is .B P .R (z) = ( P_1(z), P_2(z), P_3(z), ..., P_N_P(z) ) .br P_i(z) = ( _i(z), _i(z) , x_i(z), _xi(z), y_i(z), _yi(z) ) and .B F .R (z, .B P .R (z),a_L(x,y,z)) is the left-hand-side term of the equation of motion ()-(). With this equation, We can obtain the predictor of the particles. .B P .R ^(p)(z_n+1) = .B P .R (z_n) + .B F .R (z_n, .B P .R (z_n),a_L(x,y,z_n)) z . Now we have the predictors of the laser field and the particles. The corrector for the particles, .B P .R ^(c)(z_n+1), which is more prousible value than predictor, can be obtain by the equation .B P .R ^(c)(z_n+1) = .B P .R (z_n) + 12 [ .B F .R (z_n, .B P .R (z_n),a_L(x,y,z_n)) + .B F .R (z_n+1, .B P .R ^(p)(z_n+1),a_L^(p)(x,y,z_n+1)) ] z and .B P .R ^(c)(z_n+1) is taken as .B P .R (z_n+1). Using this corrector of the particles .B P .R ^(c)(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^(c)_j u^j C^(c)_0= s_lm(z_n) .br C^(c)_1= 16 ( 3 s_lm(z_n) - 6s_lm(z_n-1) + s_lm(z_n-2) + 2s_lm^(c)(z_n+1) ) .br C^(c)_2= 12 ( -2 s_lm(z_n) + s_lm(z_n-1) + s_lm^(c)(z_n+1) ) .br C^(c)_3= 16 ( -3 s_lm(z_n) + 3s_lm(z_n-1) - s_lm(z_n-2) + s_lm^(c)(z_n+1) ). s_lm^(c)(z) are put into the integral of the equation () and are analytically integrated, then we get the corrector,a_lm^(c)(z_n+1), which is more plausible value than the predictor a_lm^(p)(z_n+1), a_lm^(c)(z_n+1) = ( -_lm z ) a_lm(z_n) + ( _j=0^3 C^(c)_j I_lmj ) z . This a_lm^(c)(z_n+1) is taken as a_lm(z_n+1). a_L(x,y,z_n+1) can be obtaned with the inverse Fourier transformation from a_lm(z_n+1). ELFIN can get the values a_lm(z_n+1), .B P .R (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), a_lm(z_n) and .B P .R (z_n). At z=0 which is the initial point of the calculation, we do not have the history of s_lm(z) but we set can set s_lm(z)=0 for z<0. 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. .NH 1 CREATION OF ELECTRON BEAM .XS CREATION OF ELECTRON BEAM .XE .PP 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 f(x,_x,y,_y) dx d_x dy d_y = N_n (-B^n) dx d_x dy d_y, B = ( x-x x )^2 + ( _x-_x _x )^2 + ( y-y y )^2 + ( _y-_y _y )^2 , 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 /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. .NH 1 INITIAL LASER FIELD .XS INITIAL LASER FIELD .XE .PP 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) = a_L0 1+ iz-l_fz_R ( -1 1+ iz-l_fz_Rx^2+y^2w_0^2 ) .br ,where _L z_R = w_0^2, a_L0 = P2r_emc^2 1c^2 _Lw_0 .br = 4.823382413 10^-5 P[W] _Lw_0. z_R is Rayleigh length. The initial laser field on the mesh point is a_L(x_p,y_q,0). .NH 1 FREE SPACE FOR LASER FIELD .XS FREE SPACE FOR LASER FIELD .XE .PP 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 . The laser field propagating along with z axis in the free space obeys the equations ( -i2 k_L_^2 + ddz ) a_L(x,y,z) = 0 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() is the quation() with s=0. With Gaussian transformation, ( c x' .br y' .br z'-z_0' ) = ( c xz+ib .br yz+ib .br ^2 z-z_0(z_0+ib)(z+ib) ) .br g(x,y,z) = 1z+ib ( -i k_L x^2 + y^22 (z+ib) ) .br v(x',y',z') = u(x,y,z)g(x,y,z) .br here , is the complex constant and b and z_0 are the real constatnts, the equation () is transformed to ( -i2 k_L_^2' + z ) v(x',y',z') = 0 . 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 ( c x' .br y' .br z' ) = ( c xz .br yz .br z_1z ( z - z_1 ) ) .br g(x,y,z) = 1z ( -i k_Lx^2 + y^22 z ) .br v(x',y',z') = u(x,y,z)g(x,y,z) . 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'=z_1 Lz_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 (), 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 + z ) a_lm(z') = 0 and thier solutions are a_lm(z'=L) = ( -_lm L ) a_lm(z'=0) .br The laser field a_L(x,y,z'=L) can be obtained from a_lm(z'=L) by inverse Fourier transformation (). .NH 1 MIRROR REFRECTION .XS MIRROR REFRECTION .XE .PP 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) = r(x,y)^22 .br r(x,y) = (x-x)^2 + (y-y)^2 . The contribution of the tilt is z_tilt(x,y) = -x _x - y _y. 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) .br = (x-x)^2 + (y-y)^22 -x _x -y _y. The phase difference (x,y) of the laser field at point (x,y) is (x,y) = k_L d(x,y) .br = -2 k_L ( (x-x)^2 + (y-y)^22 -x _x -y _y ) ,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 C(x,y) = R(x,y) (i(x,y)) .br = R(x,y) ( 2i k_L ( -(x-x)^2+(y-y)^22 + x_x + y_y ) ) 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) where a_L^IN is the incident laser field. .NH 1 DRIFT SPACE FOR ELECTRON .XS DRIFT SPACE FOR ELECTRON .XE .PP The motion of electrons in a drift space of length L is x^out= x^in + _x^in L .br y^out= y^in + _y^in L .br where superscript in and out show that the values with them are entrance and exit of this section, respectively. .NH 1 DISPERSIVE SECTION .XS DISPERSIVE SECTION .XE .PP 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 - cv L_d ) = k_L L_d2^2 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= k_L L_d4^2 . 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+e^2m^2 c^2 L _0^L ( _0^z' B_y(z\*U) dz\*U ) dz' ] . 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 - k_L L_d2 ^2 - k_L L2( _x^in^2 + _y^in^2 ) .br x^out= x^in + _x^in L .br y^out= y^in + _y^in L .br where superscript in and out show that the values with them are entrance and exit of this section, respectively. . bibitemTang-Sprangle-85 For example, C. M. Tang and P. A. Sprangle, IEEE J. QUantum Electron. QE-21(1985),970. C. M. Tang, P. Sprangle, A. Ting and B. Hafizi , J.Appl.PHys. 66(4),15(1989) B. D. McVey, "THREE-DIMENSIONAL SIMULATIONS OF FREE ELECTRON LASER PHYSICS," Nucl. Instr. and Meth. A250(1986)449. E. A. Sziklas and A. E. Siegman, Proc. IEEE 62,410(1974). 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.