Three-Dimensional Numerical Simulations of FEL ’ s by the Transverse Mode Spectral Method

We present a numerical method for solving a three-dimensional radiation field by decomposing the radiation field into transverse modes that satisfy the free space wave equation and the boundary conditions. The formalism includes betatron oscillations in a realizable wiggler, finite emittance, and energy spread in an amplifier or oscillator configuration. This formalism can be generalized to calculate sidebands resulting from the synchrotron oscillations. The transverse mode spectral method is compared to the finite difference method, spectral method, and transform spectral method. Examples from a computer code using Gaussian-Hermite expansion are given. We show self-focusing of the radiation beam due to gain and refraction of the FEL.

5) Sideband instabilities in an oscillator can lead to pulse breakup and affect the quality of the radiation beam.
The analytical and numerical methods employed in the study of the transverse variation of the wave equation have taken the approaches of the transverse mode spectral method [1]-[5], the transform spectral method [6], [7], the spectral method [2], the finite difference method [8], [9], and the Lienard-Wiechert potential method [ 101.The different methods refer to the techniques used to evaluate the 0 : A , term in the wave equation where AR is the vector potential of the radiation field.All three kinds of spectral methods involve representing the solution of the radiation field as a truncated series of known functions of the independent variables.
The transverse mode spectral method decomposes the radiation field into a truncated series of a complete set of orthogonal functions which satisfy the free space wave equations with the appropriate boundary conditions.One takes advantage of the transverse mode properties to reduce the wave equation to a set of first-order differential equations.
For the spectral method and the transverse spectral method, the series takes on the form of a transform, and the wave equation reduces to a simple form in the transformed variables, i.e., 0 : A, has an analytical form in the transform space.The transform spectral method requires numerically taking the transform of the driving current.The spectral methods describe the current in terms of the Lagrangian variables, and the current term can be evaluated analytically.
In Section 11, we will outline the transverse mode spectral method in general, i.e., not specifying the form of the transverse modes.This formulation will include slow transverse motion and betatron oscillations of the electrons in the realistic wiggler, finite emittance, energy spread, and self-consistent axial particle dynamics.This method conserves energy.The formalism is extended to the study of sideband formation [l], [ 2 ] , [ll], [ 121 formation on a long electron pulse in the FEL oscillator.
In Section 111, we outline the spectral and the transverse spectral methods.The three different spectral methods are compared in this section.
The advantages of the transverse mode spectral method are: 1) free space wave propagation, finite size mirrors, and apertures can be handled analytically; 2) it is easy to include transverse particle motion exactly; 3 ) this method lends itself to analytical and semi-analytical solutions and U.S. Government work not protected by U.S. copyright can provide physical insight for many problems; and 4) the The orthogonality condition is transverse boundary conditions are included automatically in the waveguide mode expansion.(4) In Section IV, we apply the transverse mode spectral method to an example with self-focusing properties.In this where the integration is over the appropriate radiation example, the electron beam focuses due to gain and domain.

C
In this paper, we will consider only the linearly polar-Substituting the expression of the radiation field into the ized wiggler since that is the most common wiggler field.wave equation, and making the slowly varying approxi-The formulation for a circularly polarized wiggler requires mation, we obtain only minor modifications.The realistic magnetic field will be expressed in terms of the vector potential of the wiggler We also include a dc accelerating electric field [14] J J Jz rnoc2)A,(z).
Now we need to evaluate the current.The fluid-like electron density is E,,&, for the purpose of efficiency enhancement.For con-D ( x , y , z , p x , p y , p,, t) = do(xo, yo, P , , ~, $o) venience, we define a dimensionless parameter eDc = 1 e I ED,/rnoc2.

6(x -X) sg,jj)
The radiation field will be expressed as exp [i(kzut)] where The variables with " -" and the variable T are functions * (XP, -PJ 6(p, -P,).Gl,m = gl,,neiei~m represents the x component of the complex transverse mode (which may be zero) A&, y, z ) = [A,@, y, z)l exp (@(x, y, z ) ) , and k = u / c .We use lower the time the electrons enter the wiggler at z = 0, and T is the time that the electron reaches position z.The function do is the initial phase space density distribution, including the effects of emittance and energy spread.
The fluid-like beam density can be defined as case to denote the normalized parameters of the radiation field, such as a, (x, y, z ) = ( l e l / ~r n o c 2 > A , ( x , y , z ) and .D(X, y> z, Px, Py, pz, 9 . (9) wave equation and the appropriate boundary con-We will define the effective area of the electron beam to ditions, e., be where no = n(x, y , z , t)Imax.
The current density can be defined as The final form of the wave equation is is the ensemble average, the normalization is  and   8, is the mean axial velocity.

B. Particle Dynamics
The mean transverse particle position and momentum, in general, should be integrated along with the axial particle momentum equations.Under certain circumstances, analytical forms for the transverse particle trajectories are sufficient.
Taking the betatron oscillations as an example, the particle motion in the x direction, without any additional external focusing, can be written as The particle motion in the y direction exhibits betatron oscillations.Taking k,jj << 1, the particle dynamics in The particle motion in the z direction is best written in terms of the equation for the total relativistic gamma factor and the equation governing the phase where is the degree of taper for the efficiency enhancement schemes, includes the effect resulting from the combination of betatron oscillation and the contoured B,. field, and Equations (12), (16), and (17) form the self-consistent set of equations for the radiation field and the particle dynamics.

C. Conservation of Energy
We can show that this formulation conserves energy.TO do this, we rewrite (16) 41 where x = h x / w ( z ) , = &y/w(z), H! is the Hermite polynomial of the Ith order, j-= ( zz,)/zo, z, is the axial location of the minimum waist, w(z) = wo( 1 + j-2)''2, w o is the spot size, and zo = &k/2 is the Rayleigh length.
2) When the electron beam and the resonator system are axially symmetrical, and when the betatron oscillations are not included in the model of the linearly polarized wiggler, then the radiation can be assumed to be axially symmetrical.The Gaussian-Laguerre modes are the choice for the expansion, and each term is expressed as 1 = G0,,S, where "2 and L i is the Laguerre polynomial 3) In a low gain oscillator, the best choice would be the eigenmodes of the resonator cavity [18] including the effects of apertures, finite size mirrors, and resonator losses.The eigenmodes in turn can be written in terms of the appropriate Gaussian-Laguerre or Gaussian-Hermite modes.The advantage of going to the resonator cavity modes is that the losses for the higher order resonator cavity modes increases rapidly with the cavity mode number, such that only a few modes need to be kept in the calculation.In addition, there is no need to calculate the radiation outside the FEL wiggler.

4)
Closed waveguides can be used to confine the radiation beam over a long distance compared to a Rayleigh length.If the waveguide is not rectangular, the vector potentials of the waveguide modes generally have both x and y components.For the analysis presented here, we have to assume that the transverse guide dimensions are large compared to a radiation wavelength.The theory is also restricted to low-order, low-loss modes whose propagation constants k + are nearly equal to the plane-wave value, i.e., Re@/,,) << k .For wavelengths that are comparable to the waveguide diameters or' when multimode analysis is desired, this method has to be modified slightly to include the different frequencies and axial wavenumbers associated with the different waveguide modes.The formulation should be similar to that described in Section 11-F with sidebands.

E. Computation Speed
The cost of computation per axial step is spent mainly for the radiation field.The number of operations is a,N, ( L + l)(M + 1) to integrate (12), where N , is the total number of electrons and a' is a numerical constant.The particle orbit equations require a2(L + 1)(M + 1) N,/No number of operations to evaluate the radiation field a&, g , z ) and a 3 N e number of operations for the integration where No is the number of electrons per ponderomotive wave and a2 and a3 are numerical constants.In general, ( 16) is faster than (20), which requires a4(L + 1) (M + 1)N, number of operations.

F. Frequency Sidebands ,and Pulse Slippage
The electron oscillation in the ponderomotive potential well is called the synchrotron oscillation.We define Kso to be the synchrotron wavenumber of the electron traveling exacly along the z axis: When the intensity of the radiation becomes high, such that the electrons make half a bounce in the wiggler, then sideband frequencies can grow.One-dimensional simulations [ 121 and three-dimensional simulations [ 11, [2] show that the amplitude of the radiation field becomes chaotic and the quality of the radiation field becomes degraded.
In order to understand the effect of sidebands on wavefront curvatures, the three-dimensional formulation is necessary.We assume periodic boundary conditions for the radiation field of length l e b , which is chosen to be many times the amplitude modulation distance and pulse slippage distance.If the length of the electron beam is much longer than le,,, only a section of a long electron beam needs to be modeled.If the electron beam is shorter than le,,, then the whole pulse shape must be included.
The radiation field can be written as The principle behind the spectral and transform spectral methods is similar to the transverse mode spectral method.The major difference comes from the representation of the radiation field, for example, the Fourier series in Cartesian geometry and the Hankel transform in axially cylindrical geometry.
Let us take the Fourier series for the illustration.The domain for the radiation field i s -0 J 2 5 x I D x / 2 and -0 J 2 I y 5 0 , / 2 and the boundary conditions are periodic.The slowly varying amplitude of the radiation field is written as The method described above is straightforward.Since the number of electrons that must be used to model the electron beam with finite emittance and finite length is large, any general numerical method, including this method, requires a considerable amount of computation time.For limited computational resources, it is necessary to form simplified models.This formulation lends itself towards semi-analytical models of finite length electron pulses [l], [2], [18] where some of the integrals in the expression for ensemble average can be evaluated analytically.

(33)
When the current J, is formulated in Eulerian variables, for the x and y coordinates then the right-hand side of (33) has to be evaluated numerically, and the method is called the transform spectral method.If the current is formulated in Lagrangian variables (xO, yo,  px,O, P , , ~, pz,o, $, , ) then ( 3 3 ) can be reduced to The particle dynamics equations (16), (17) are also applicable here.For a Fourier series expansion, one can show that the energy is conserved for the spectral method.
In general, the Eulerian formulation of the current is inferior to the Lagrangian formulation.The Eulerian formulation is not convenient to use for the study of betatron effects because the numerical transforms of (33) require the current to be evaluated at prespecified grids.This problem can be reduced if the number of grids used across the electron beam is large, which in turn requires a larger number of terms in the expansion.For the Hankel trans- form, there is the additional problem of numerical errors resulting from numerical integration with a finite number of grids in an infinite integration domain.
The estimate for the speed of computation for the spectral method has the same form as in Section 11-E.The values for the coefficients a l , a 2 , and a3 depend on the expansion, and are smaller for the Fourier transform and larger for the Hankel transform.
The transverse mode spectral method, in some instances, is better than both the spectral and transform spectral methods because of the following properties.The propagation of the radiation field through apertures and the reflection and transmission of the radiaton field at mirrors can be evaluated in terms of matrix multiplication of the amplitudes of the normal modes.The spectral and transform spectral methods require an additional step in converting the radiation field into Gaussian modes.Finally, the transverse mode spectral method is easier to implement for calculations in complicated waveguide geometries because the boundary conditions are automatically included.

IV. SELF-FOCUSING EXAMPLE AND DISCUSSION
A three-dimensional code employing Gaussian-Hermite expansion of the radiation field is applied to a linearly polarized wiggler.The following examples will show selffocusing properties of the FEL's.
The electron beam is assumed to have a Gaussian profile and a radius of 2.25 X cm.The current is 50 A and the energy is 109.6 MeV.The wiggler has a magnetic field B , = 6.3 kG, a period of 2.4 cm, and a length L, = 6 m.The incident 0.5 pm radiation at the entrance of the wiggler is a Gaussian TEMoo mode with a minimim waist wo = 0.06 cm located at z = L,/2 and a power of 5.  150, 300, 450, and 600 cm plotted from -4wo to 4wo in both the x and y plane.In this example, the self-focusing phenomenon is not only due to refraction, but also gain.The peak amplitude in Fig. 2(a)-(e) ,are 4.5 X 1.1 X lop2, 4.6 X 4.1 X lop2, and 2.6 X respectively.At 150 cm, the laser beam is already narrowed down significantly.The corresponding phase front is shown in Fig. 3(a).The focusing due to the electron beam is manifested by the small mound at the center.
The beam radius remains at roughly the same size from 150 cm to the end of the wiggler.Due to electrons oscillating away from the bottom of the ponderomotive potential well, the radiation loses energy after 4 m.In [13] and [ 141, it was shown analytically that FEL radiation is not only governed by gain and diffraction, but also by refraction when operating in the trapped particle mode.For this case, the resonant phase is at the origin, i.e., sin GR = 0.Even though the radiation is losing energy after z = 4 m, the radiation beam around the electrons still focuses, while the radiation further away from the electrons defocuses, as shown in Fig. 2(e).This is in qualitative agreement with the theory given by [13] and [14].Fig. 3(b) shows the focusing properties of the wavefront.It is important that the radiation beam does not defocus when the radiation is losing energy due to the bouncing of the electrons in the ponderomotive potential well.
The FEL radiation does not always undergo focusing, for example, the electron beam could dig a hole in the laser amplitude profile when the FEL is operating with loss in the low gain regime and large frequency mismatch.. Next, we give an example where the center of the electron beam does not travel down the center of the wiggler due to irregularities in the magnetic wiggler field.The center of the electron beam is assumed to execute a slow sinusoidal motion ycntr = reb sin ( ~z l L , ) .The radiation guiding properties of the e beam are still clearly evident (Fig. 4).The peak amplitude is 2.4 X lo-'.The motion of the center of the electron beam, however, causes significant distortion of the radiation beam.
In summary, this paper outlines one method of solving the three-dimensional FEL radiation field self-consistently with the electron dynamics.The advantages of this method are: 1) the boundaries in the transverse directions are included automatically; 2) it is easy to incorporate transverse particle motion; 3) free space prapagation, finite size mirrors, and apertures can be handled analytically; and 4) this method lends itself to analytical and semi-analytical solutions.
The disadvantage of this method is the increase in the computation time if a large number of electrons are used in the radial direction, and if a large number of modes are required.
The numerical examples illustrate the self-focusing property of the FEL.Under appropriate conditions, the laser radiation maintains a roughly constant radius radiation beam.We also showed that the irregularities in the wiggler field can cause significant distortion of the radia-  tion beam.The implications of this for the FEL design are important.For electron beams with good emittance, it is possible to focus the beam to a small area and operate the FEL in the high gain regime, and not have to worry about diffraction.In this case, it is important to design the resonator cavity including the FEL gain model.

( 8 )
of the Lagrangian varaibles XO, yo, p , , ~, p y , 0 , p , , ~, $0, z where xo, yo are the initial transverse positions at the en-L M trance of the wiggler z = 0, P , , ~, py,o are the initial trans-= C C ~[ , , ( z ) ii,m(x, y, z ) eiP'~~~~z, verse momentums at z = 0, pZ,o is the initial axial momentum spread at z = 0, 1L0 = -uto is the initial phase / = 0 m=O &l,m = G l > m i x + F l , m z y of the electron in the ponderomotive potential well, to is is the general form of the transverse mode, Al,m = IAl,,Ieic'*" is the complex amplitude of the normal modes, .
in terms of the amplitude of each mode: * I a1, ,(z) I e sin -CDC, (20) -Im(P/.dzwhere +l,m = $ + + C P ~, , ~ is the phase associated with each mode.We have assumed here that Re(Pl,,>z << 2 a .The spatial rate of change of the total electron beam energy is (moc2dT/dz) = -DFl(z) A,.(z) .1Al,&) / e -w / .d z 1,m .( g / , m ( -f , 9,;' sin $l,m ) -(moc2eDc>, (21) where D = (o,,/8a)(wi/c2)k.The spatial rate of change in the radiation energy is Applying (12) to (22), we can show that the energy is conserved, i.e., the energy gained by the radiation is equal to the energy lost by the electrons plus the energy supplied by the axial dc electric field: D. The Choice of the Transverse Modes of the Wave Equation The choice of the normal modes, however, will depend on the geometry of the FEL of interest.A few examples of the possible choices of the normal eigenmodes are given below.1) Gaussian-Hermite expansion allows the most flexibility in the modeling of the FEL physics of interest in an open resonator.The expressions for varying function of position following the radiation pulse, k = (1 + /3,0)(y$(l + K 2 ) ) k, is the resonant wavenumber, Ak,, = (n -1)6k, 6k = 27r/leb, and 4 = zct.A(x, y , 4, t ) is a slowly varying function of position and time following the radiation pulse.Following the same procedure outlined in Sections II-A and II-B, we obtain the self-consistent set of equations for al,m,n, and $.The independent Lagrangian variables are xo3 Y O , pX,o, P ~, o , P ~, o , $0, * o , t , and average, \EO = (2.ir/leb)to, and 0 I E,, I leb are the axial location of the electron in the electron pulse at t = 0, $ = jb(k,(ct') + k)b, dt ' -ckt + $o, 2 = j b 6, dt' + to, = z" -ct = \kO/6k -s, and s = ct(1 + K 2 ) / 2 y i is the pulse slippage distance.

6 AxialFig. 1 .
Fig. 1.A plot of the spatial power gain as a function of axial distance.

8 X
In this example, the electron beam radius is much smaller than the spot size, i.e., reb << wo.If the radiation has the spot size of the electron beam, its Rayleigh length would be only 32 cm.A plot of the power gain is shown in Fig.1.The FEL is operating in the high gain regime, and the radiation saturates at 4 m.io5 w.

Fig. 4 .
Fig. 4. Plot of the amplitude of the radiation field in the n and y plane for Addison-Wesley, 1980, pp.147-174.the electron beam that deviated from the axis, The tick marks are sepa-[I21J.c.Goldstein and w.B. ColSon, ''Control Of optical Pulse modurated by a distance wo.lation due to the sideband instability in free electron lasers," in Proc.Int.Con$ Lasers'82, R .C. Powell, Ed.McLean, VA: STS Press,