MITHRA 1.0: A full-wave simulation tool for Free Electron Lasers
MMITHRA 1.0
A F
ULL - WAVE S IMULATION T OOL for Free Electron Lasers
Arya Fallahi, Alireza Yahaghi and Franz X. K ¨artner
Group of Ultrafast Optics and X-ray SourcesDESY-Center for Free Electron Laser Science (CFEL) a r X i v : . [ phy s i c s . c o m p - ph ] D ec Mithra, also spelled Mithras, Sanskrit Mitra, in ancient Indo-Iranian mythology, is the god of light, whosecult spread from India in the east to as far west as Spain, Great Britain, and Germany. The first writtenmention of the Vedic Mitra dates to 1400 BC. His worship spread to Persia and, after the defeat of the Persiansby Alexander the Great, throughout the Hellenic world. In the 3rd and 4th centuries AD, the cult of Mithra,carried and supported by the soldiers of the Roman Empire, was the chief rival to the newly developing religionof Christianity. The Roman emperors Commodus and Julian were initiates of Mithraism, and in 307 Diocletianconsecrated a temple on the Danube River to Mithra, “Protector of the Empire.”According to myth, Mithra was born, bearing a torch and armed with a knife, beside a sacred stream andunder a sacred tree, a child of the earth itself. He soon rode, and later killed, the life-giving cosmic bull, whoseblood fertilizes all vegetation. Mithra’s slaying of the bull was a popular subject of Hellenic art and becamethe prototype for a bull-slaying ritual of fertility in the Mithraic cult.As god of light, Mithra was associated with the Greek sun god, Helios, and the Roman Sol Invictus. He is oftenpaired with Anahita, goddess of the fertilizing waters. Source: Encyclopaedia Britannica ontents
MESH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.2
BUNCH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.3
FIELD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.4
UNDULATOR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.5
FEL-OUTPUT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 ist of Tables ist of Figures E y at 110 µ m distance from the bunch center and (b) the total radiated power calculated at110 µ m distance from the bunch center in terms of the traveled undulator length. . . . . . . . . . . . . . . . . . . . . 384.2 Snapshots of the radiated field profile taken at (a) x = 0 and (b) z = 60 µm plane and (c) the bunch profile viewedfrom the x axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.3 Convergence study for the different involved parameters in the considered FEL simulation: (a) n , (b) ∆ t b , (c) l z , (d) l x = l y , (e) ∆ z and (f) ∆ x = ∆ y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394.4 The total radiated power calculated at 110 µ m distance from the bunch center in terms of the traveled undulator length(a) with and without space-charge consideration and (b) various lengths of the bunch with space-charge assumption. . 404.5 (a) The total radiated power measured at 80 µ m distance from the bunch center in terms of the traveled undulatorlength and (b) the bunch profile at 12 m from the undulator begin. . . . . . . . . . . . . . . . . . . . . . . . . . . . 434.6 The total radiated power calculated at 110 µ m distance from the bunch center in terms of the traveled undulatorlength compared for two cases of an optical and static undulator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.7 (a) Electric field of the generated radiation in front of the bunch, (b) the total radiated power measured at 82 nmdistance from the bunch center in terms of the traveled distance, (c) the same radiation power for various harmonicorders, and (d) bunching factor of the considered bunch during the ICS interaction. . . . . . . . . . . . . . . . . . . 484.8 The total radiated power measured at 82 nm distance from the bunch center in terms of the traveled distance for animaginary bunch where each electron is represented by a cloud of 1000 particles. . . . . . . . . . . . . . . . . . . . 494 hapter 1 Introduction
Free Electron Lasers (FELs) are currently serving as promising and viable solutions for the generation of radiation in the wholeelectromagnetic spectrum ranging from microwaves to hard X-rays [1, 2, 3]. Particularly, in portions of spectrum where commonsolutions like lasers and other electronic sources do not offer efficient schemes, FEL based devices attract considerable attention andinterest. For example, soft and hard X-ray radiation sources as well as THz frequency range are parts of the spectrum where FELsources are widely used. In the optical regime, lasers currently serve as the most popular sources, where radiation is generated andamplified based on the stimulated emission. More accurately, the excited electrons of the gain medium emit coherent photons whenchanging the energy level to the ground state [4]. Since the energy bands of different gain media are fixed curves determined by thematerial atomic lattice, there are only specific wavelengths obtainable from lasers operating based on stimulated emission in a gainmedium. In contrast, there exist vacuum electronic devices like gyrotrons, klystrons and travelling wave tubes (TWT), in which freeelectrons travelling along a certain trajectory transform kinetic energy to an electromagnetic wave [5]. Although, these sources areusually not as efficient as medium based lasers, their broadband operations make them promising in portions of the spectrum whereno gain media is available.In a free electron laser, relativistic electrons provided from linear accelerators travel through a static undulator and experience awiggling motion. The undulator performance is categorized into two main regimes: ( i ) in a short undulator, each electron radiates asan independent moving charge, which yields an incoherent radiation of electron bunch. Therefore, the radiation power and intensityis linearly proportional to the number of electrons. ( ii ) For long interaction lengths, the radiated electromagnetic wave interacts withthe bunch and the well-known micro-bunching phenomenon takes place. Micro-bunching leads to a periodic modulation of chargedensity inside the bunch with the periodicity equal to the radiation wavelength. This effect results in a coherent radiation scalingwith the square of the bunch numbers. Coherent X-ray have shown unprecedented promises in enabling biologists, chemists andmaterial scientists to study various evolutions and interactions with nanometer and sub-nanometer resolutions [6].Owing to the desire of hard X-ray FEL machines for electrons with ultrarelativistic energies (0.5-1 GeV), these sources areusually giant research facilities with high operation costs and energy consumption. Therefore, it is crucial and additionally veryuseful to develop sophisticated simulation tools, which are able to capture the important features in a FEL radiation process. Suchtools will be very helpful for designing and optimizing a complete FEL facility and additionally useful for detailed investigation ofimportant effects. The last decade had witnessed extensive research efforts aiming to develop such simulation tools. As a result,various softwares like Genesis 1.3 [7], MEDUSA [8], TDA3D [9, 10], GINGER [11], PERSEO [12], EURA [13], RON [14], FAST[15], PlaRes [16] and Puffin [17] are developed and introduced to the community. However, all the currently existing simulationsoftwares are usually written to tackle special cases and therefore particular assumptions or approximations have been consideredin their development [18]. Some of the common approximations in FEL simulation are tabulated in Table 1.1. The main goal inthe presented software is the analysis of the FEL interaction without considering any of the above approximation. The outcome ofthe research and effort will be a sophisticated software with heavy computation loads. Nonetheless, it provides a tool for testing thevalidity of various approximations in different operation regimes and also a reliable approach for preparing the final design of a FELfacility.Besides the wide investigations and studies on the conventional X-ray FELs, recently research efforts have been devoted tobuilding compact X-ray FELs, where novel schemes for generating X-ray radiations in a so-called table-top setup are examined andassessed. Various research topics such as laser-plasma wake-field acceleration (LPWA) [19, 20, 21], laser plasma accelerators (LPA)[22, 23], laser dielectric acceleration (LDA) [24] and THz acceleration [25, 26], pursue the development of compact acceleratorscapable of delivering the desired electron bunches to FEL undulators. Besides such attempts, one promising approach to make acompact undulator is using optical undulators, where the oscillations in an electromagnetic wave realize the wiggling motion of5 HAPTER 1. INTRODUCTION 6
Table 1.1: Common approximations in modelling free electron laser radiationcode name approximation1D FEL steady state wiggler-average slow wave forward no space-charge slicetheory approximation electron motion approximation waveGENESIS 1.3 — (cid:88) (optional) (cid:88) (cid:88) (cid:88) — (cid:88) MEDUSA — (cid:88) (optional) — (cid:88) (cid:88) — (cid:88) TDA3D — (cid:88) (cid:88) (cid:88) (cid:88) — no time-domainGINGER — — (cid:88) (cid:88) (cid:88) — —PERSEO — — — — (cid:88) (cid:88) —EURA — — (cid:88) (cid:88) (cid:88) — —FAST — — (cid:88) (cid:88) — — (cid:88)
Puffin — — — — (cid:88) (cid:88) —the electrons [27]. Many of the approximations in Table 1.1, which sound reasonable for static undulators are not applicable forstudying an optical undulator radiation. In this regime, due to the various involved length-scales and remarkable impact of theparameter tolerances, having access to a rigorous and robust FEL simulation tool is essential.One of the difficulties in the X-ray FEL simulation stems from the involvement of dramatically multidimensional electromagneticeffects. Some of the nominal numbers in a typical FEL simulation are: • Size of the bunch: ∼ fs or µ m • Undulator period: ∼ cm • Undulator length: ∼ − m • Radiation wavelength: ∼ − nmComparing the typical undulator lengths with radiation wavelengths immediately communicates the extremely large space for thevalues. This in turn predicts very high computation costs to resolve all the physical phenomena, which is not practical even with theexisting supercomputer technology. In order to overcome this problem, we exploit Lorentz boosted coordinate system and implementFinite Difference Time Domain (FDTD) [28] method combined with Particle in Cell (PIC) simulation in the electron rest frame. Thiscoordinate transformation makes the bunch size and optical wavelengths longer and shortens the undulator period. Interestingly,these very different length scales transform to values with the same order after the coordinate transformation. Consequently, thelength of the computation domain is reduced to slightly more than the bunch length making the full-wave simulation numericallyfeasible. We comment that the simulation of particle interaction with an electromagnetic wave in a Lorentz boosted framework isnot a new concept. The advantage of this technique is already demonstrated and widely used in the simulation of plasma-wakefieldacceleration [29, 30, 31]. However, to the best of our knowledge, it has never been used to simulate a FEL mechanism, which iswhat we are aiming in this study.The presented manual shows how one can numerically simulate a complete FEL interaction using merely Maxwell equations,equation of motion for a charged particle, and the relativity principles. In chapter 2, the whole computational aspects of the soft-ware, including the Finite Difference Time Domain (FDTD), Particle In Cell (PIC), current deposition, Lorentz boosting, quantityinitialization, and parallelization, are described in detail. The implementation is explained in a way suitable for a graduate studentto start writing the code on his own. Chapter 3 provides a reference card for a software user to get familiar with MITHRA and therequired parameters for performing the simulations. Finally, in chapter 4, different examples of free electron lasers are analyzedand the results are presented in conjunction with some discussions. As a new software entering the FEL community, we aim tokeep updating this material with new implementations and examples. In this regard, any assistance and help from the users of thissoftware will be highly appreciated. hapter 2 Methodology
In this chapter, we present the detailed formalism of Finite Difference Time Domain - Particle In Cell (FDTD/PIC) method in theLorentz boosted coordinate system. There are many small still very important considerations in order to obtain reliable results, whichconverge to the real values. For example, the method for electron bunch generation, particle pusher algorithm and computationalmesh truncation need particular attention.
FDTD is perhaps the first choice coming to mind for solving partial differential equations governing the dynamics of a system.Despite its simple formulation and second order accuracy, there are certain features in this method like explicit time update and zeroDC fields, which makes this method a superior choice compared to other algorithms [28]. FDTD samples the field in space and timeat discrete sampling points and represents the partial derivatives with their discrete counterparts. Subsequently, update equations arederived based on the governing differential equation. Using these updating equations, a time marching algorithm is acquired whichevaluates the unknown functions in the whole computational domain throughout the simulation time. In the following, we start withthe Helmholtz equation which is the governing partial differential equation for our electromagnetic problem.
The physics of electromagnetic wave and its interaction with charged particles in free space is mathematically formulated throughthe well-known Maxwell’s equations: ∇ × E = − ∂ B ∂t (2.1) ∇ × B = µ J + µ ε ∂ E ∂t (2.2) ∇ · E = − ρε (2.3) ∇ · B = 0 (2.4)These equations in conjunction with the electric current equation J = ρ v ( v is the charge velocity) and the Lorentz force equation: F = q ( E + v × B ) (2.5)are sufficient to describe wave-electron interaction in free space. Moving free electrons introduce electric current which enters intothe Maxwell’s equations as the source. Electric and magnetic fields derived from these equations are subsequently employed inthe Lorentz force equation to determine the forces on the electrons, which in turn determine their motions. As it is evident fromthe above equations, there are two unknown vectors ( E and B ) to be evaluated, meaning that six unknown components should beextracted from the equations. However, since these two vectors are interrelated and specially because there is no magnetic monopolein the nature ( ∇ · B = 0 ), one can recast Maxwell’s equations in a Helmholtz equation for the magnetic vector potential ( A ) and aHelmholtz equation for the scalar electric potential ( ϕ ): ∇ A − c ∂ ∂t A = − µ J (2.6)7 HAPTER 2. METHODOLOGY 8 ∇ ϕ − c ∂ ϕ∂t = − ρε , (2.7)where c = 1 / √ µ ε is the light velocity in vacuum. In the derivation of above equations, the Lorentz gauge ∇ · A = − c ∂ϕ∂t isused. The original E and B vectors can be obtained from A and ϕ as: B = ∇ × A (2.8) E = − ∂ A ∂t − ∇ ϕ (2.9)In addition to the above equations, the charge conservation law written as ∇ · J + ∂ρ∂t = 0 , (2.10)should not be violated in the employed computational algorithm. This is the main motivation for seeking proper current depositionalgorithms in the FDTD/PIC methods used for plasma simulations. It is immediately observed that the equations (2.6), (2.7), (2.10)and the Lorentz gauge introduce an overdetermined system of equations. In other words, once a current deposition is implementedthat automatically satisfies the charge conservation law, the Lorentz gauge will also hold, provided that the scalar electric potential( ϕ ) is obtained from (2.7). However, due to the space-time discretization and the interpolation of quantities to the grids, a suitablealgorithm that holds the charge conservation without violating energy and momentum conservation does not exist. The approachthat we follow in MITHRA is using the discretized form of (2.6) and (2.7) with the currents and charges of electrons (i.e. macro-particles) as the source and solving for the vector and scalar potential. It can be shown that if current density ( J ) and charge density( ρ ) are interpolated similarly to the vertices of each cell, a proper discretization of current density based on positions of the macro-particles holds the charge conservation criterion. To obtain the fields E and B at the grid points, we use the momentum conservinginterpolation, which will be explained in the upcoming sections. In cartesian coordinates, a vector Helmholtz equation is written in form of three uncoupled scalar Helmholtz equations. Therefore,it is sufficient to apply our discretization scheme only on a typical scalar Helmholtz equation: ∇ ψ − c ∂ ψ∂t = ζ , where ψ standsfor A l ( l ∈ { x, y, z } ); and ζ represents the term − µ J l . Let us begin with the central-difference discretization scheme for variouspartial differential terms of the scalar Helmholtz equation at the point ( i ∆ x, j ∆ y, k ∆ z, n ∆ t ) . In the following equations, ψ ni,j,k denotes the value of the quantity ψ at the point ( i ∆ x, j ∆ y, k ∆ z ) and time n ∆ t . The derivatives are written as follows: ∂ ∂x ψ ( x, y, z, t ) (cid:39) ψ ni +1 ,j,k − ψ ni,j,k + ψ ni − ,j,k (∆ x ) (2.11) ∂ ∂y ψ ( x, y, z, t ) (cid:39) ψ ni,j +1 ,k − ψ ni,j,k + ψ ni,j − ,k (∆ y ) (2.12) ∂ ∂z ψ ( x, y, z, t ) (cid:39) ψ ni,j,k +1 − ψ ni,j,k + ψ ni,j,k − (∆ z ) (2.13) ∂ ∂t ψ ( x, y, z, t ) (cid:39) ψ n +1 i,j,k − ψ ni,j,k + ψ n − i,j,k (∆ t ) . (2.14)Combining these four equations, one obtains the value of ψ at instant ( n + 1)∆ t in terms of its value at n ∆ t and ( n − t : ψ n +1 i,j,k = − ψ n − i,j,k + α ψ ni,j,k + α ψ ni +1 ,j,k + α ψ ni − ,j,k + α ψ ni,j +1 ,k + α ψ ni,j − ,k + α ψ ni,j,k +1 + α ψ ni,j,k − + α ζ ni,j,k where the coefficients α , . . . , α are obtained from: α = 2 (cid:34) − (cid:18) c ∆ t ∆ x (cid:19) − (cid:18) c ∆ t ∆ y (cid:19) − (cid:18) c ∆ t ∆ z (cid:19) (cid:35) , α = ( c ∆ t ) ,α = α = (cid:18) c ∆ t ∆ x (cid:19) , α = α = (cid:18) c ∆ t ∆ y (cid:19) , α = α = (cid:18) c ∆ t ∆ z (cid:19) (2.15) HAPTER 2. METHODOLOGY 9
The term ζ ni,j,k is the magnitude of the source term at the time n ∆ t , which is calculated from the particle motions. Usually, one needsa finer temporal discretization for updating the equation of motion compared to electromagnetic field equations. If the equation ofmotion is discretized and updated with ∆ t b = ∆ t/N time steps, the term ζ ni,j,k will be written in terms of the value after each N update: ζ ni,j,k = − µ J l ( n ∆ t ) = − µ ρ ( n ∆ t ) r n +1 / − r n − / ∆ t . (2.16)As observed in the above equation, the position of particles are sampled at each n + 1 / time step, which later should be consideredfor updating the scalar potential. This assumption also results in the calculation of charge density at n +1 / time steps, which shouldbe averaged for obtaining ρ ( n ∆ t ) . It is well-known that the FDTD formulation for discretizing the Helmholtz equation suffers from the so-called numerical dispersion.More accurately, the applied discretization leads to the wave propagation with a speed different from (lower than) the vacuum speedof light. This may impact the FEL simulation results particularly during the saturation regime, owing to the important role played bythe relative phase of electrons with respect to the radiated light. Therefore, careful scrutiny of this effect and minimizing its impactis essential for the goal pursued by MITHRA.To derive the equation governing such a dispersion, we assume a plane wave function for ψ ( x, y, z, t ) = e − j ( k x x + k y y + k z z − ωt ) in the discretized Helmholtz equation. After some mathematical operations, the following equation is obtained for the dispersionproperties of central-difference scheme: sin ( k x ∆ x/ x ) + sin ( k y ∆ y/ y ) + sin ( k z ∆ z/ z ) = sin ( ω ∆ t/ c ∆ t ) . (2.17)This equation is evidently different from the vacuum dispersion relation, which reads as k x + k y + k z = ω c . (2.18)Comparison of the two equations shows that the dispersion characteristics are similar, if and only if ∆ x → , ∆ y → , ∆ z → ,and ∆ t → . Another output of the dispersion equation is the stability condition, which is referred to as Courant-Friedrichs-Lewy(CFL) condition [28]. The spatial and temporal discretization should be related such that the term ω obtained from equation (2.17)has no imaginary part, i.e. sin ( ω ∆ t/ < . This implies that c ∆ t < (cid:113) sin ( k x ∆ x/ x ) + sin ( k y ∆ y/ y ) + sin ( k z ∆ z/ z ) . (2.19)The right hand side of the above equation has its minimum when all the sinus functions are equal to one, which leads to the stabilitycondition for the central-difference scheme: ∆ t < c (cid:113) x ) + y ) + z ) . (2.20)As mentioned above, for the FEL simulation, it is very important to maintain the vacuum speed of light along the z direction(Throughout this document z is the electron beam and undulator direction). More accurately, if k x = k y = 0 , k z = ω/c should bethe solution of the dispersion equation. However, this solution is obtained if and only if ∆ t = ∆ z/c , which violates the stabilitycondition. To resolve this problem, various techniques are developed in the context of compensation of numerical dispersion. Here,we take advantage from the non-standard finite difference (NSFD) scheme to impose the speed of light propagation along z direction[32, 33].The trick is to consider a weighted average along z for the derivatives with respect to x and y , which is formulated as follows: ∂ ∂x ψ ( x, y, z, t ) (cid:39) ¯ ψ ni +1 ,j,k − ψ ni,j,k + ¯ ψ ni − ,j,k (∆ x ) (2.21) ∂ ∂y ψ ( x, y, z, t ) (cid:39) ¯ ψ ni,j +1 ,k − ψ ni,j,k + ¯ ψ ni,j − ,k (∆ y ) , (2.22) HAPTER 2. METHODOLOGY 10 with ¯ ψ ni,j,k = A ψ ni,j,k − + (1 − A ) ψ ni,j,k + A ψ ni,j,k +1 . (2.23)Such a finite difference scheme leads to the following dispersion equation: (cid:0) − A sin ( k z ∆ z/ (cid:1) (cid:18) sin ( k x ∆ x/ x ) + sin ( k y ∆ y/ y ) (cid:19) + sin ( k z ∆ z/ z ) = sin ( ω ∆ t/ c ∆ t ) . (2.24)It can be shown that if the NSFD coefficient A is larger than 0.25, and (cid:112) (∆ z/ ∆ x ) + (∆ z/ ∆ y ) < , a real ω satisfies the abovedispersion equation for ∆ t = ∆ z/c . This time step additionally yields k z = ω/c , for k x = k y = 0 .The value we chose for A in MITHRA is obtained from A = 0 . (cid:18) . z/ ∆ x ) + (∆ z/ ∆ y ) (cid:19) . (2.25)The update equation can then be written as ψ n +1 i,j,k = − ψ n − i,j,k + α (cid:48) ψ ni,j,k + α (cid:48) ( A ψ ni +1 ,j,k − + (1 − A ) ψ ni +1 ,j,k + A ψ ni +1 ,j,k +1 ) + α (cid:48) ( A ψ ni − ,j,k − + (1 − A ) ψ ni − ,j,k + A ψ ni − ,j,k +1 )+ α (cid:48) ( A ψ ni,j +1 ,k − + (1 − A ) ψ ni,j +1 ,k + A ψ ni,j +1 ,k +1 ) + α (cid:48) ( A ψ ni,j − ,k − + (1 − A ) ψ ni,j − ,k + A ψ ni,j − ,k +1 )+ α (cid:48) ψ ni,j,k +1 + α (cid:48) ψ ni,j,k − + α (cid:48) ζ ni,j,k . (2.26)where the coefficients α (cid:48) , . . . , α (cid:48) are obtained from: α (cid:48) = 2 (cid:34) − (1 − A ) (cid:18) c ∆ t ∆ x (cid:19) − (1 − A ) (cid:18) c ∆ t ∆ y (cid:19) − (cid:18) c ∆ t ∆ z (cid:19) (cid:35) , α (cid:48) = ( c ∆ t ) ,α (cid:48) = α (cid:48) = (cid:18) c ∆ t ∆ x (cid:19) , α (cid:48) = α (cid:48) = (cid:18) c ∆ t ∆ y (cid:19) , α (cid:48) = α (cid:48) = (cid:18) c ∆ t ∆ z (cid:19) − A (cid:18) c ∆ t ∆ x (cid:19) − A (cid:18) c ∆ t ∆ x (cid:19) . (2.27)To guarantee a dispersion-less propagation along z direction with the speed of light the update time step is automatically calculatedfrom the given longitudinal discretization ( ∆ z ), according to ∆ t = ∆ z/c . Usually, due to high energy of particles in a FEL process, the FEL simulations neglect the space-charge effects by considering ϕ (cid:39) [16]. However, this is an approximation which we try to avoid in MITHRA. To account for space-charge forces, one needs to solvethe Hemholtz equation for scalar potential, i.e. (2.7). For this purpose, the same formulation as used for the vector potential isutilized to update the scalar potential. Nonetheless, since the position of particles are sampled at t + ∆ t/ instants, the obtainedvalue for ϕ n corresponds to the scalar potential at ( n + 1 / t . This point should be particularly taken into consideration, whenelectromagnetic fields E and B are evaluated. In order to simulate the FEL problem, we consider a cube as our simulation domain. The absorbing boundary condition is alsoconsidered for updating the scalar electric potential ϕ at the boundaries. Therefore, we introduce the parameter ξ , which denoteseither ψ or ϕ . The six boundaries of the cube are supposed to be at: x = ± l x / , y = ± l y / and z = ± l z / . In the following,the process for implementing Mur absorbing boundary conditions (ABCs) of the first and second order in MITHRA are discussed.We only present the formulation for the boundary conditions at z = ± l z / . The process to extract the equations for the other fourboundaries will be exactly similar. First Order ABCs:
The partial differential equations implying first order ABCs at z = ± l z / are: ∓ ∂ ξ∂z∂t − c ∂ ξ∂t = 0 (2.28)The discretized version for different terms appearing in the above equation reads as: HAPTER 2. METHODOLOGY 11 • At z = − l z / k = 0) ∂ ξ∂z∂t (cid:39) t (cid:32) ξ n +1 i,j, − ξ n +1 i,j, ∆ z − ξ n − i,j, − ξ n − i,j, ∆ z (cid:33) (2.29) c ∂ ξ∂t (cid:39) c (cid:32) ξ n +1 i,j, − ξ ni,j, + ξ n − i,j, ∆ t + ξ n +1 i,j, − ξ ni,j, + ξ n − i,j, ∆ t (cid:33) (2.30) • At z = + l z / k = K = l z / ∆ z ) ∂ ξ∂z∂t (cid:39) t (cid:32) ξ n +1 i,j,K − ξ n +1 i,j,K − ∆ z − ξ n − i,j,K − ξ n − i,j,K − ∆ z (cid:33) (2.31) c ∂ ξ∂t (cid:39) c (cid:32) ξ n +1 i,j,K − ξ ni,j,K + ξ n − i,j,K ∆ t + ξ n +1 i,j,K − − ξ ni,j,K − + ξ n − i,j,K − ∆ t (cid:33) (2.32)Combining these equations, one obtains the boundary value of ξ at instant ( n + 1)∆ t in terms of its values at n ∆ t and ( n − t : • At z = − l z / k = 0) ξ n +1 i,j, = β ξ n − i,j, + β ξ ni,j, + β ξ n − i,j, + β ξ ni,j, + β ξ n +1 i,j, (2.33) • At z = + l z / k = K = l z / ∆ z ) ξ n +1 i,j,K = β ξ n − i,j,K + β ξ ni,j,K + β ξ n − i,j,K − + β ξ ni,j,K − + β ξ n +1 i,j,K − (2.34)where: β = β = c ∆ t − ∆ zc ∆ t + ∆ z , β = β = 2∆ zc ∆ t + ∆ z , β = − (2.35) Second Order ABCs:
The partial differential equations implying second order ABCs at z = ± l z / are: ∓ ∂ ξ∂z∂t − c ∂ ξ∂t − c ∂ ξ∂x − c ∂ ξ∂y = 0 (2.36)The discretized version for different terms appearing in the above equation reads as: • At z = − l z / k = 0) ∂ ξ∂z∂t (cid:39) t (cid:32) ξ n +1 i,j, − ξ n +1 i,j, ∆ z − ξ n − i,j, − ξ n − i,j, ∆ z (cid:33) (2.37) c ∂ ξ∂t (cid:39) c (cid:32) ξ n +1 i,j, − ξ ni,j, + ξ n − i,j, ∆ t + ξ n +1 i,j, − ξ ni,j, + ξ n − i,j, ∆ t (cid:33) (2.38) c ∂ ξ∂x (cid:39) c (cid:18) ξ ni +1 ,j, − ξ ni,j, + ξ ni − ,j, ∆ x + ξ ni +1 ,j, − ξ ni,j, + ξ ni − ,j, ∆ x (cid:19) (2.39) c ∂ ξ∂y (cid:39) c (cid:18) ξ ni,j +1 , − ξ ni,j, + ξ ni,j − , ∆ y + ξ ni,j +1 , − ξ ni,j, + ξ ni,j − , ∆ y (cid:19) (2.40) HAPTER 2. METHODOLOGY 12 • At z = + l z / k = K = l z / ∆ z ) ∂ ξ∂z∂t (cid:39) t (cid:32) ξ n +1 i,j,K − ξ n +1 i,j,K − ∆ z − ξ n − i,j,K − ξ n − i,j,K − ∆ z (cid:33) (2.41) c ∂ ξ∂t (cid:39) c (cid:32) ξ n +1 i,j,K − ξ ni,j,K + ξ n − i,j,K ∆ t + ξ n +1 i,j,K − − ξ ni,j,K − + ξ n − i,j,K − ∆ t (cid:33) (2.42) c ∂ ξ∂x (cid:39) c (cid:18) ξ ni +1 ,j,K − ξ ni,j,K + ξ ni − ,j,K ∆ x + ξ ni +1 ,j,K − − ξ ni,j,K − + ξ ni − ,j,K − ∆ x (cid:19) (2.43) c ∂ ξ∂y (cid:39) c (cid:18) ξ ni,j +1 ,K − ξ ni,j,K + ξ ni,j − ,K ∆ y + ξ ni,j +1 ,K − − ξ ni,j,K − + ξ ni,j − ,K − ∆ y (cid:19) (2.44)Combining these equations, one obtains the boundary value of ξ at instant ( n + 1)∆ t in terms of its values at n ∆ t and ( n − t : • At z = − l z / k = 0) ξ n +1 i,j, = γ ξ n − i,j, + γ ξ ni,j, + γ ξ n − i,j, + γ ξ ni,j, + γ ξ n +1 i,j, + (2.45) γ ξ ni +1 ,j, + γ ξ ni − ,j, + γ ξ ni,j +1 , + γ ξ ni,j − , + γ ξ ni +1 ,j, + γ ξ ni − ,j, + γ ξ ni,j +1 , + γ ξ ni,j − , • At z = + l z / k = K = l z / ∆ z ) ξ n +1 i,j,K = γ ξ n − i,j,K + γ ξ ni,j,K + γ ξ n − i,j,K − + γ ξ ni,j,K − + γ ξ n +1 i,j,K − + (2.46) γ ξ ni +1 ,j,K − + γ ξ ni − ,j,K − + γ ξ ni,j +1 ,K − + γ ξ ni,j − ,K − + γ ξ ni +1 ,j,K + γ ξ ni − ,j,K + γ ξ ni,j +1 ,K + γ ξ ni,j − ,K where: γ = γ = c ∆ t − ∆ zc ∆ t + ∆ z , γ = γ = ∆ z (cid:0) − ( c ∆ t/ ∆ y ) − ( c ∆ t/ ∆ x ) (cid:1) c ∆ t + ∆ z , γ = − γ = γ = γ = γ = ( c ∆ t/ ∆ x ) ∆ z c ∆ t + ∆ z ) , γ = γ = γ = γ = ( c ∆ t/ ∆ y ) ∆ z c ∆ t + ∆ z ) (2.47)Particular attention should be devoted to the implementation of Mur second order absorbing boundary condition at edges andcorners. Separate usage of the above equations for second order case encounters problems in the formulation. On one hand, unknownvalues at grid points outside the computational domain appears in the equations, and on the other a system of overdeterminedequations will be obtained. The solution to this problem is to discretize all the involved boundary conditions at the center of the cubes(for corners) or squares (for edges). A simple addition of the obtained equations cancels out the values outside the computationaldomain and returns the desired value meeting the considered absorbing boundary condition. Particle in cell (PIC) method is the standard algorithm to solve for the motion of particles within an electromagnetic field distribution.The method takes the time domain data of the fields E and B and updates the particle position and momentum according tothe Lorentz force equation (2.5). We comment that the electromagnetic fields in the motion equation are the total fields in thecomputational domain, which in a FEL problem is equivalent to the superposition of undulator field, radiated field and the seededfield in case of a seeded FEL problem. Often considering all the individual particles involved in the problem ( ∼ − particles)leads to high computation costs and long simulation times. The clever solution to this problem is the macro-particle assumption,through which an ensemble of particles ( ∼ − particles) are treated as one single entity with charge to mass ratio equal to theparticles of interest, which are here electrons. The relativistic equation of motion for electron macro-particles then reads as ∂∂t ( γm v ) = − e ( E + v × B ) , and ∂ r ∂t = v , (2.48)where r and v are the position and velocity vectors of the electron, e is the electron charge and m is its rest mass. γ stands for theLortenz factor of the moving particle. HAPTER 2. METHODOLOGY 13
Δx ΔzΔy zxy(iΔx,jΔy,kΔz)δx δz δy
Figure 2.1: Schematic illustration of the parameters used to locate a particle within the computational domain.
There are numerous update algorithms proposed for the time domain solution of (2.48), including various Runge-Kutta and finitedifference algorithms. Among these methods, Boris scheme has garnered specific attention owing to its interesting peculiarity whichis being simplectic. Simplectic update algorithms are update procedures which maintain the conservation of any parameter in theequation which obey a physical conservation law. Since in a FEL problem effect of the magnetic field on a particle motion plays themost important role, using a simplectic algorithm is essential to obtain reliable results. This was the main motivation to choose theBoris scheme for updating the particle motion in MITHRA.We sample the particle position at times m ∆ t b , which is represented by r m and the particle normalized momentum at times ( m − )∆ t b which is written as γ β m − / . Then, by having r m and γ β m − / as the known parameters and E mt and B mt as the total field values imposed on the particle at instant m ∆ t , the values r m +1 and γ β m +1 / are obtained as follows: t = γ β m − / − e ∆ t b E mt mc t = t + α t × B mt t = t + t × α B mt α B mt · B mt γ β m +1 / = t − e ∆ t b E mt mc r m +1 = r l + c ∆ t b γ β m +1 / (cid:113) γ β m +1 / · γ β m +1 / (2.49)with α = − e ∆ t b / (2 m √ t · t ) . E mt = E mext + E m and B mt = B mext + B m are total fields imposed on the particle, whichare equal to the superposition of the radiated field with the external fields, i.e. the undulator or the seed fields. In order to figure outthe derivation of the equations (2.49), the reader is referred to [34, 35]. As seen from the above equations, the electric and magneticfields at time m ∆ t b and the position r of the particle are needed to update the motion. In the next section, the equations to extractthese values from the computed values of the magnetic and scalar potential are presented. Note that to achieve a certain precisionlevel, the required time step in updating the bunch properties ( ∆ t b ) is usually much smaller than the time step for field update ( ∆ t ).In MITHRA, there exists the possibility for setting different time steps for PIC and FDTD algorithms. As described in section 2.1, the propagating fields in the computational domain are evaluated by solving the Helmholtz equationfor the magnetic vector potential, i.e. (2.6). To update the particle position and momentum, one needs to obtain the field values E m and B m from the potentials A and ϕ . For this purpose, the equations (2.8) and (2.9) need to be discretized in a consistentmanner to provide the accelerating field with lowest amount of dispersion and instability error. First, the values of magnetic andscalar potentials at t + ∆ t/ are used to evaluate the electromagnetic fields at the cell vertices. Subsequently, the field values are HAPTER 2. METHODOLOGY 14 interpolated to the particle location for updating the equation of motion. An important consideration at this stage is compatibleinterpolation of fields from the cell vertices with the interpolations used for current and charge densities. Similar interpolationalgorithms should be followed to cancel the effect of self-forces on particle motion.Using the equation (2.8), the magnetic field B ni,j,k at cell vertex ( i, j, k ) is calculated as follows: B xni,j,k = 12 (cid:32) A zni,j +1 ,k − A zni,j − ,k y − A yni,j,k +1 − A yni,j,k − z + A zn +1 i,j +1 ,k − A zn +1 i,j − ,k y − A yn +1 i,j,k +1 − A yn +1 i,j,k − z (cid:33) , (2.50) B yni,j,k = 12 (cid:32) A xni,j,k +1 − A xni,j,k − z − A zni +1 ,j,k − A zni − ,j,k x + A xn +1 i,j,k +1 − A xn +1 i,j,k − z − A zn +1 i +1 ,j,k − A zn +1 i − ,j,k x (cid:33) , (2.51) B zni,j,k = 12 (cid:32) A yni +1 ,j,k − A yni − ,j,k x − A xni,j +1 ,k − A xni,j − ,k y + A yn +1 i +1 ,j,k − A xn +1 i − ,j,k x − A xn +1 i,j +1 ,k − A xn +1 i,j − ,k y (cid:33) . (2.52)Similarly, equation (2.9) is employed to evaluate the electric field at the cell vertices. The electric field E ni,j,k is obtained from thefollowing equations: E xni,j,k = (cid:32) − A xn +1 i,j,k − A xni,j,k ∆ t − ϕ ni +1 ,j,k − ϕ ni − ,j,k x (cid:33) , (2.53) E yni,j,k = (cid:32) − A yn +1 i,j,k − A yni,j,k ∆ t − ϕ ni,j +1 ,k − ϕ ni,j − ,k y (cid:33) , (2.54) E zni,j,k = (cid:32) − A zn +1 i,j,k − A zni,j,k ∆ t − ϕ ni,j,k +1 − ϕ ni,j,k − z (cid:33) . (2.55)Suppose that a particle resides at the cell ijk with the grid point indices shown in Fig. 2.1. As illustrated in Fig.2.1, thedistance to the corner ( i ∆ x, j ∆ y, k ∆ z ) is assumed to be ( δx, δy, δz ) . We use a linear interpolation of the fields from the ver-tices to the particle position to calculate the imposed field. If ς denotes for a component of the electric or magnetic field, i.e. ς ∈ { E x , E y , E z , B x , B y , B z } , one can write ς p = (cid:88) I,J,K (cid:18)
12 + ( − I (cid:12)(cid:12)(cid:12)(cid:12) − δx ∆ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) (cid:18)
12 + ( − J (cid:12)(cid:12)(cid:12)(cid:12) − δy ∆ y (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) (cid:18)
12 + ( − K (cid:12)(cid:12)(cid:12)(cid:12) − δz ∆ z (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ς i + I,j + J,k + K , (2.56)where I , J , and K are equal to either 0 or 1, producing the eight indices corresponding to the eight corners of the mesh cell. Once the position and momentum of all the particles over the time interval ∆ t is known, one needs to couple the pertinent currentsinto the Helmholtz equation (2.6). As described before, this coupling over time is implemented through the equation (2.16). Theremaining question is how to evaluate the related currents on the grid points, i.e. the method for performing an spatial interpolation.To maintain consistency, we should use a similar interpolation scheme as used for the field evaluation. This assumption leads to thefollowing equation for spatial interpolation. ρ pi + I,j + J,k + K = ρ (cid:18)
12 + ( − I (cid:12)(cid:12)(cid:12)(cid:12) − δx ∆ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) (cid:18)
12 + ( − J (cid:12)(cid:12)(cid:12)(cid:12) − δy ∆ y (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) (cid:18)
12 + ( − K (cid:12)(cid:12)(cid:12)(cid:12) − δz ∆ z (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) (2.57)where ρ is the charge density attributed to each macro-particle, namely q/ (∆ x ∆ y ∆ z ) . ρ pi,j,k is the charge density at the grid point ( i, j, k ) due to the moving particle p in the computational mesh cell (Fig. 2.1a). I , J , and K are equal to either 0 or 1, which producethe eight indices corresponding to the eight corners of the mesh cell. The total charge density ρ i,j,k will be a superposition of all thecharge densities due to the moving particles of the bunch. We have removed the superscripts corresponding to the time instant, toavoid the confusion due to different time marching steps ∆ t and ∆ t b . The above interpolation is carried out at each update step ofthe field values. One can consider the above interpolation equations as a rooftop charge distribution centered at the particle positionand expanding in the regions ( − ∆ x < x < ∆ x, − ∆ y < y < ∆ y, − ∆ z < z < ∆ z ) . Eventually, the equation (2.16) is used tocalculate the corresponding current densities. HAPTER 2. METHODOLOGY 15
Laboratory coordinate system Bunch rest frame:
Figure 2.2: Schematic illustration of the Lorentz boosting to transform the problem from the laboratory frame to the bunch restframe.
The previous two sections on FDTD and PIC algorithms present a suitable and efficient framework for the computation of interactionbetween charged particles and propagating waves. However, the initial conditions are always required for a complete determinationof the problem of interest. For a FEL simulation, the initial conditions corresponding to the FEL input are given to the FDTD/PICsolver. For example, in case of a SASE (Self Amplified Spontaneous Emission) FEL, the initial fields are zero and there is noexcitation entering the computational domain, whereas for a seeded FEL, an outside excitation should be considered entering thecomputational domain. The explanation of how such initializations are implemented in MITHRA is the goal in this section.One novel feature of the method, followed here, is the solution of Maxwell’s equations in the bunch rest frame. It can be shownthat a proper coordinate transformation yields the matching of all the major parameters in a FEL simulation, namely bunch length,undulator period, undulator length, and radiation wavelength. Fig. 2.2 schematically describes the advantage of moving into thebunch rest frame. In a typical FEL problem, the FEL parameter ρ F EL is about − . Therefore, simulation of FEL interactionwith a bunch equal to the cooperation length of the FEL ( L c = λ l / (4 πρ F EL ) , with λ l being the radiation wavelength) requires asimulation domain only 100 times larger than the wavelength. This becomes completely possible with the today computer technologyand constitutes the main goal of MITHRA. In this section, the main basis for Lorentz boosting the simulation coordinate is describedfirst. Afterwards, the relations for evaluating the undulator fields in the Lorentz boosted framework are presented. Finally, theelectron bunch initialization in the Lorentz-boosted framework is discussed. It is known from the FEL thoery that a bunch with central Lorentz factor equal to γ moves in an undulator with an average Lorentzfactor equal to γ = γ/ (cid:112) (1 + K / , where K = eBλ u / (2 πmc ) is the undulator parameter determining the amplitude of thewiggling motion. Consequently, a frame moving with normalized velocity β = (cid:112) − /γ is indeed the bunch rest frame, wherethe volume of the computational domain stays minimal. Transforming into this coordinate system necessitates tailoring the bunchand undulator properties. For this purpose, the Lorentz length contraction, time dilation and relativistic velocity addition need to beemployed.In MITHRA, the input parameters are all taken in the laboratory frame and the required Lorentz transformations are carried outbased on the bunch energy. The required transformations for the computational mesh are as the following: ∆ z = ∆ z (cid:48) γ , (2.58) ∆ t = ∆ t (cid:48) /γ , (2.59) ∆ t b = ∆ t (cid:48) b /γ , (2.60)where the prime sign stands for the quantities in the laboratory frame. The quantities without prime are values in the bunch restframe, which are used in the FDTD/PIC simulation. With the consideration of the above transformations, the length of the totalcomputational domain along the undulator period and the total simulation time is also transformed similarly.In addition to the data for the computational mesh, the properties of the electron bunch also changes after the Lorentz boosting.This certainly affects the bunch initialization process which is thoroughly explained in the next section. An electron bunch inMITHRA is initialized and characterized by the following parameters:(i) Mean electron position: (¯ x b , ¯ y b , ¯ z b ) , HAPTER 2. METHODOLOGY 16
Laboratory coordinate system z'x' y'
Figure 2.3: Schematic illustration of the undulator in the lab frame and the definition of the coordinates.(ii) Mean electron normalized momentum: ( γβ x , γβ y , γβ z ) ,(iii) RMS value of the electron position distribution: ( σ x , σ y , σ z ) ,(iv) RMS value of the electron normalized momentum distribution: ( σ γβ x , σ γβ y , σ γβ z ) .As mentioned previously, the above parameters are entered by the user in the laboratory frame. To transform the given values to thebunch rest frame the position related parameters are changed as ¯ x b = ¯ x (cid:48) b , ¯ y b = ¯ y (cid:48) b , ¯ z b = ¯ z (cid:48) b γ (1 − ¯ β (cid:48) z β ) , σ x = σ (cid:48) x , σ y = σ (cid:48) y , σ z = σ (cid:48) z γ (1 − ¯ β (cid:48) z β ) . (2.61)To transfer the momentum related quantities, we assume that the main contribution to the Lorentz factor is the momentum along z direction or the undulator period. In other words, ( γβ x , γβ y , γβ z ) = γ ( ¯ β x , ¯ β y , ¯ β z ) , with γ = 1 / (cid:112) − ¯ β z . Similarly, theRMS values can also be written as ( σ γβ x , σ γβ y , σ γβ z ) = γ ( σ β x , σ β y , σ β z ) . Using the relativistic velocity transformation [36], thetransformation equations for the above values are found as follows: γ = γ (cid:48) γ (1 − ¯ β (cid:48) z β ) , (2.62) ( ¯ β x , ¯ β y , ¯ β z ) = ( ¯ β (cid:48) x , ¯ β (cid:48) y , (cid:112) − /γ ) , (2.63) ( σ γβ x , σ γβ y , σ γβ z ) = ( σ (cid:48) γβ x , σ (cid:48) γβ y , σ (cid:48) γβ z ) γ (1 − ¯ β (cid:48) z β ) . (2.64)Equations (2.58)-(2.64) provide a sufficient set of equations to perform the Lorentz boost to the bunch rest frame. The utilized FDTD/PIC algorithm solves the Maxwell’s equation coupled with the motion equation of an ensemble of particles.Therefore, in addition to the field values, particle initial conditions should also be initialized. For a SASE FEL problem, the initialfield profile is zero everywhere, whereas for a seeded FEL the initial seed should enter the computational domain through theboundaries. In both cases, the external field which is the undulator field should separately be initialized.
Undulator Field:
By solving the Laplace equation for the magnetic field, the undulator field in the laboratory frame is found to be as the following(Fig. 2.3) [1]: B (cid:48) x = 0 ,B (cid:48) y = B cosh( k u y (cid:48) ) sin( k u z (cid:48) ) , (2.65) B (cid:48) z = B sinh( k u y (cid:48) ) cos( k u z (cid:48) ) , where B is the maximum transverse field of the undulator. To calculate the undulator field in the bunch rest frame, first the positionis transformed to laboratory frame ( x (cid:48) , y (cid:48) , z (cid:48) ) through the Lorentz boost equations. Afterwards, the field is evaluated using the HAPTER 2. METHODOLOGY 17 equation (2.65). Eventually, these fields are transformed back into the bunch rest frame. Following this procedure, the undulatorfield in the bunch rest frame is obtained as B x = 0 , E x = 0 B y = γ B cosh( k u y ) sin( k u l z ) , E y = 0 (2.66) B z = B sinh( k u y ) cos( k u l z ) , E z = γ β cB cosh( k u y ) sin( k u l z ) , with l z = γ ( z (cid:48) + β ct − z ) and z being the initial distance between the bunch and undulator in the bunch rest frame.An important consideration in the initialization of undulator field is the entrance region of the undulator. A direct usage of theequation (2.66) with zero field for z (cid:48) < causes an abrupt variation in the particles motion, which results in a spurious coherentradiation. In fact, in a real undulator, there exists fringing fields at the undulator entrance, which remove any abrupt transitionin the undulator field and consequently the particle radiations [37]. To the best of our knowledge, the fringing fields are alwaysmodeled numerically and there exists no analytical solution for the problem. Here, we approximate the fringing fields by a graduallydecreasing magnetic field in form of a Neumann function. The coefficients in the function are set such that the particles do not gainany net transverse momentum and stay in the computational domain as presumed. The undulator field for z (cid:48) < in the laboratoryframe is obtained as the following: B (cid:48) x = 0 ,B (cid:48) y = B cosh( k u y (cid:48) ) k u z (cid:48) e − ( k u z (cid:48) ) / , (2.67) B (cid:48) z = B sinh( k u y (cid:48) ) e − ( k u z (cid:48) ) / , The same transformations as in (2.66) can be used to approximate the fringing field values in the bunch rest frame.
Seed Field:
External excitation of free electron laser process using a seed mechanism has proved to be advantageous in terms of output spec-trum, photon flux and the required undulator length [1, 38]. Such benefits has propelled the proposal of seeded FEL schemes. Tosimulate such a mechanism, MITHRA uses the TF/SF (total-field/scattered-field) technique to introduce an external excitation intothe computational domain. When seeding is enabled by having a non-zero seed amplitude, the second and third points (after theboundary points) constitute the scattered and total field boundaries, respectively. Therefore, during the time marching process, aftereach update according to equation (2.26) the excitation terms are added to the fields at TF/SF boundaries. For example for the TF/SFboundaries close to z = z min plane, the field values to be used in the next time steps are obtained as the following:SF boundary: ψ (cid:48) n +1 i,j,k = ψ n +1 i,j,k + A ( α (cid:48) f ni +1 ,j,k +1 + α (cid:48) f ni − ,j,k +1 + α (cid:48) f ni,j +1 ,k +1 + α (cid:48) f ni,j − ,k +1 ) + α (cid:48) f ni,j,k +1 , TF boundary: ψ (cid:48) n +1 i,j,k = ψ n +1 i,j,k − A ( α (cid:48) f ni +1 ,j,k − + α (cid:48) f ni − ,j,k − + α (cid:48) f ni,j +1 ,k − + α (cid:48) f ni,j − ,k − ) − α (cid:48) f ni,j,k − , (2.68)where f ni,j,k is the excitation value at time n ∆ t and position ( i ∆ x, j ∆ y, k ∆ z ) . The excitation value is calculated based on theimposed seed fields, which are usually either a plane wave or a Gaussian beam radiation. Position and momentum initialization:
As described previously, the evolution of the electron bunch is always simulated by following the macro-particle approach, wherean ensemble of particles are represented by one sample particle. This typically reduces the amount of computation cost for updatingthe bunch properties by three or four orders of magnitude. Due to the high sensitivity of a FEL problem to the initial conditions,correct and proper initialization of these macro-particles play a critical role in obtaining reliable results. In computational accel-erator physics, different approaches are introduced and developed for bunch generation. Some examples are random generation ofparticles, mirroring macro-particles at different phases to prevent initial average bunching factors, and independent random fillingof different coordinates to prevent unrealistic correlations [39]. Among all the different methods, using the sophisticated methodsto load the bunch in a ”quasi-random” manner seem to be the most appropriate solutions. The Halton or Hammersley sequences,as generalizations of the bit-reverse techniques, are implemented in MITHRA for particle generation. These sequences comparedto random based filling of the phase space avoid the appearance of local clusters in the bunch distribution. In addition, the uniformfilling of the phase space prevents initial bunching factor of the generated electron bunch, making it well-suited for FEL simulations.
HAPTER 2. METHODOLOGY 18
For details on the nature of Halton sequences, the reader is referred to the specialized documents. Here, we only present theimplemented algorithm to generate the required sequence of numbers filling the interval [0 , . The following C++ code is integratedinto MITHRA which produces N < uncorrelated sequences of M numbers in the interval [0 , : std::vector
12 ) + y m . (2.69)Such a distribution is used when a bunch with uniform current profile ( z distribution of particles) is to be initialized. On the otherhand, a 1D Gaussian distribution is needed when radiation of a bunch with Gaussian current profile is modelled. To generate buncheswith Gaussian distribution, we employ Box-muller’s theory to extract a sequence of numbers with Gaussian distribution from twouncorrelated uniform distributions. Based on this theory, a 1D Gaussian distribution y with average y m and deviation width y s isfound by the following transformation: y = y s (cid:112) − x cos(2 πx ) + y m . (2.70)Similar to the undulator fields, an abrupt variation in the bunch profile results in an unrealistic coherent scattering emission(CSE), which happens if the uniform bunch distribution is directly initialized from equation (2.69). CSE is avoided by imposingsmooth variations in the particle distribution. For this purpose, we follow the procedure proposed in [40] and [39]. A small Gaussianbunch with the same density as the real bunch and a width equal to an undulator wavelength is produced. The lower half of the bunch(particles with smaller z ) is transferred to the tail and the other half is placed at the head of the uniform bunch. Hence, a uniformcurrent profile with smooth variations at its head and tail is created.The transverse coordinates of the bunches are initialized using 2D distributions. In MITHRA, a 2D Gaussian distribution isassumed for transverse coordinates. To generate such a distribution, two independent sets of numbers x and x are generated HAPTER 2. METHODOLOGY 19 node node node node node node
MEMORY
BIOS node node node node node node memory memory memory memory memory memory
BIOS ( a ) ( b ) Figure 2.4: Schematic illustration of the (a) shared and (b) distributed memory parallelization schemes.based on Halton sequence. The desired 2D Gaussian distribution with average position ( y m , y m ) and total deviation ( y s , y s ) isproduced as the following: y = y s (cid:112) − x cos(2 πx ) + y m , and y = y s (cid:112) − x sin(2 πx ) + y m . (2.71)Such algorithms are similarly used to generate the distribution in particle momenta. The only difference is that for initializinga distribution in momentum merely Gaussian profiles are considered in transverse and longitudinal coordinates. The method tointroduce these bunch types are described in the next chapter. Bunching factor:
Free electron laser radiation should start from a nonzero initial radiation. This radiation can be in form of an initial seed field, initialmodulation in the bunch, or the radiation from bunch shot noise. In the current version of MITHRA, we have implemented thefirst two types. The implementation of shot noise is postponed to the next versions of the software. The implementation of seedingthrough an external excitation using TF/SF boundaries was described in 2.3.2. Here, we explain how an initial bunching factor, < e jk u z > , is introduced to the electron bunch profile.For this purpose, a small variation δz is applied to a particle distribution, which is generated using the above formulations. δz for each particle is obtained from δz = ξγ k u b f sin(2 ξγ k u z ) , (2.72)where b f is the given bunching factor of the distribution, and ξ = 1 + ¯ β z /β accounts for the change in the bunch longitudinalvelocity after entering the undulator. The introduced variation to the bunch coordinates, i.e. z → z + δz , yields a bunch with all thegiven particle and momentum distributions and the desired bunching factor, b f . The large and demanding computation cost needed for the simulation of the FEL process even in the Lorentz boosted coordinateframe necessitates solving the problem on multiple processors to achieve reasonable computation times. Therefore, efficient paral-lelization techniques should be implemented in the FDTD/PIC algorithm to develop an efficient software. Traditionally, there aretwo widely used techniques to run a computation in parallel on several processors: (1) shared memory, and (2) distributed mem-ory parallelization. In the shared memory parallelization or the so-called multi-threading technique, several processors run a codeusing the variables saved in one shared memory (Fig. 2.4a). This technique is very suitable for PIC algorithms because it avoidsthe additional costs of communicating the particle position and momenta between the processors. On the other hand, distributedmemory technique distributes the involved variables among several processors, solves the problem in each processor independentlyand communicates the required variables whenever they are called (Fig. 2.4b). The distributed memory technique is very suitable forFDTD algorithm due to the ease of problem decomposition beyond various machines. The advantage is fast reading and writing ofthe data and the possibility to share the computational load between different machines.Choosing a suitable parallelization scheme for the hybrid FDTD/PIC algorithm depends on both problem size and machine im-plementations. Therefore, both techniques are implemented in the software and finding an efficient distribution of the computationaltask among the available processors is offered to the user. More accurately, the user has the possibility to determine how manyprocessors are threaded (shared memory) and how far the multi-threaded algorithm is parallelized among separate processors. Themulti-threading parallelization in MITHRA is implemented using OpenMP package. The implementation simply requires distribut-ing the loops among different threads along with certain considerations to avoid race conditions. The total computational domain is
HAPTER 2. METHODOLOGY 20 z ... ... processor set i processor set i+1 communication region Figure 2.5: Schematic illustration of the domain decomposition used for distributed memory parallelization in MITHRAthen decomposed to several separate regions, each of them solved by one set of multi-threaded processors. These sets of processorscommunicate the required variables based on the technique visualized in Fig. 2.5.To parallelize the computation among N sets of multi-threaded processors, the whole computational domain is divided into N domains along z (undulator period) axis. In each time update of the field, the field values at the boundaries of each domain arecommunicated with the corresponding processor set. To parallelize the PIC solver, we define a communication domain which asshown in Fig. 2.5, is the region between the boundaries of each processor. After each update of the particles position, it is checkedif the particle has entered a communication domain. In case of residing in the communication region, the master processor, whichis the processor containing the particle in the previous time step, communicates the new coordinates to the slave processor, which isthe processor sharing the communication region with the master one. Through this simple algorithm, both parallelization schemesfunction simultaneously to achieve the fastest computation feasible and compatible with an available computing machine. hapter 3 User Interface
This chapter, as apparent from its name, is considered as a reference card for the MITHRA code. We aim to present the functionsand variables which can be delivered to the MITHRA software and can be handled for a FEL simulation problem. In what followsin this chapter, we introduce the defined language of MITHRA to write a compatible job file. This chapter can also be consideredas a reference for the current capabilities of MITHRA and with time will be updated with the further improvement of the softwarecapabilities.
Iron Rule: parameters that are used for the solution of a specific electromagnetic problem are delivered to the code at only onesingle location, the job file . This is indeed the only thing that the solver takes as an input parameter.It should be noted that all the parameters in the job file is given in the laboratory frame. The Lorentz boost into the bunch rest framewill be done by the software automatically.To run a job file using the software MITHRA, the following command should be written in the linux command line: • mpirun -np ”number of distributed processors” ”MITHRA object file name” ”job file name”The transferred job file to the solver contains five main sections, each one defining an essential part of the electromagneticproblem. These sections include:1. MESH:
The parameters of the FDTD solver like the computational domain, cell sizes and time steps are set in this section.2.
BUNCH:
The required data to initialize the electron bunch in the computational domain is set in this section. In addition, thedesired type of recording the bunch evolution is entered in this section by the user.3.
FIELD:
This section fulfills the same task as the previous section for the electromagnetic fields. The field initialization incase of a seeded FEL and the desired output type for the field evolution is given in this section to the software.4.
UNDULATOR:
This section introduced the different parameters of the undulator.5.
FEL-OUTPUT:
The desired data concerning the FEL radiation and how to record this data is set in this section.In the next subsections, we explain each part and the supported parameters, respectively. To write comments in your job file use thesign ”
MESH
As mentioned above, this part is dedicated to the determination of the FDTD/PIC parameters. In Fig 3.1, a typical computationdomain assumed in MITHRA is depicted. The mesh and update parameters of the solver are defined through the following elevenparameters: • length-scale is the scaling of the length and all the spatial parameters in the job file. The capability to play with lengthscales is crucial to avoid working with very large or very small numbers. • time-scale is the scaling of the time and all the temporal parameters in the job file. Similar to above, through this capabilityworking with very large or very small numbers is avoided.21 HAPTER 3. USER INTERFACE 22 zy Δz Δyl y l z xy Δxl x ( x ,y ,z ) Figure 3.1: The definition of the spatial mesh parameters in MITHRA • mesh-lengths is a three dimensional vector equal to the lengths ( l x , l y , l z ) of the computational domain (3.1) along thethree Cartesian axes. • mesh-resolution defines the length of one single grid cell or in other words the spatial discretization resolution of theFDTD mesh in the laboratory coordinate system (∆ x (cid:48) , ∆ y (cid:48) , ∆ z (cid:48) ) . • mesh-center is the position of the central point of the computational rectangle, i.e. ( x (cid:48) , y (cid:48) , z (cid:48) ) in Fig. 3.1. • total-time is the total computation time in the scale given by the time scale. This is indeed the time it takes for the electronbunch to travel through the considered undulator length. • bunch-time-step is the time step for updating the macro-particles’ coordinates in the PIC solver. • bunch-time-start is the time point when the bunch of macro-particles emerges in the computational domain. • mesh-truncation-order is the truncation order of the absorbing boundary condition in the computational domain. Thisparameter can be either 1 or 2, representing the first order and second order absorbing boundary condition. • number-of-threads is the level of multi-threading in the solver. This corresponds to the number of processors using ashared memory in one processor set. • space-charge is a boolean flag determining if the space-charge effect should be considered or not. If this flag is false,the scalar potential φ is zero throughout the calculation. Otherwise, the scalar potential is calculated using the correspondingHelmholtz equation.The format of the MESH group is:
MESH{ length-scale = < real value | METER | DECIMETER | CENTIMETER | MILLIMETER| MICROMETER | NANOMETER | ANGSTROM >time-scale = < real value | SECOND | MILLISECOND | MICROSECOND |NANOSECOND | PICOSECOND | FEMTOSECOND | ATTOSECOND >mesh-lengths = < ( real value, real value, real value) >mesh-resolution = < ( real value, real value, real value) >mesh-center = < ( real value, real value, real value) >total-time = < real value >bunch-time-step = < real value >bunch-time-start = < real value >mesh-truncation-order = < 1 | 2 >number-of-threads = < integer value >space-charge = < true | false >}
An example of the computational mesh definition looks as the following:
HAPTER 3. USER INTERFACE 23
MESH{ length-scale = MICROMETERtime-scale = PICOSECONDmesh-lengths = ( 3200, 3200.0, 280.0)mesh-resolution = ( 50.0, 50.0, 0.1)mesh-center = ( 0.0, 0.0, 0.0)total-time = 30000bunch-time-step = 1.6bunch-time-start = 0.0number-of-threads = 2mesh-truncation-order = 2space-charge = false}
Note that there are some conditions, which should be fulfilled for the numerical integrator to obtain reliable dispersion-lessresults. The software checks for these conditions before starting to solve the problem, if the conditions are violated the closest valueto the given number meeting the violated conditions will be used. Regarding the above parameters. the software checks for thestability condition (cid:112) (∆ z/ ∆ x ) + (∆ z/ ∆ y ) < , adapts the values of ∆ x and ∆ y accordingly, and finally sets the time step forfield update equal to ∆ z/c . In addition, the bunch update time step should be an integer fraction of the field time step to avoidredundant dispersion in the calculated values. Therefore, the closest value to the given bunch time step, which satisfies the abovecriterion, will be chosen. BUNCH
The section
BUNCH is the main part of the job file to establish the required data for the bunch input and output framework. Thissection consists of four groups: (1) bunch-initialization , (2) bunch-sampling , (3) bunch-visualization , and(4) bunch-profile . As apparent from the name the first group determines the set of parameters to initialize the bunch and theother three groups are dedicated to reporting the bunch evolution in different formats. In what follows, the parameters in each groupare introduced:1. bunch-initialization:
This group mainly determines the parameters whose values are needed for initializing a bunchof electrons with different types. If several bunches are present in a simulation, this group should simply be repeated in the
BUNCH section. The set of values accepted in this group include: • type is the type of the bunch to be initialized in the computational domain. There are five bunch types supported byMITHRA:(a) manual initializes charges at the points specified by the position vector. At each appearance of this type ofbunch only one single macro-particle will be initialized. Therefore, to have multiple manual initialization, the bunch-initialization group should be repeated. Using the file type is a better solution for high numberof manual inputs.(b) ellipsoid initializes charges with a(c) initializes multiple bunches on the points of a 3D crystal centered at the coordinate specified by theposition vector and extends over the space by the number vector and the considered lattice constant. Each singlebunch has a ellipsoid Gaussian property with the values read from the deviation parameters.(d) file reads a list of 6D position and momentum coordinates from a file and initializes the macro-particles corre-spondingly in the solver. • distribution determines if the initialized particle distribution should have a uniform or Gaussian current profile. • number-of-particles is the total number of particles (or macro-particles) considered in the bunch. • charge is the total charge of the bunch in one electron charge unit. • gamma is the initial mean Lorentz factor of the . • beta is the initial mean normalized velocity of the particles, if it is not determined here the value will be calculatedfrom the gamma parameter, otherwise the same beta will be used. HAPTER 3. USER INTERFACE 24 • direction is the average momentum direction of the bunch, i.e. ( β x , β y , β z ) /β . • position is the central position of the bunch. • sigma-position is the deviation in position of the bunch, i.e. ( σ x , σ y , σ z ) for Gaussian distributions. σ z half thebunch length for the uniform distribution. • sigma-momentum is the deviation in energy of the bunch, i.e. ( σ γβ x , σ γβ y , σ γβ z ) . • numbers is a parameter read only when the bunch type is a type. It is the number of bunch replicationin the three directions. • lattice-constants is a parameter read only when the bunch type is a type. It is the length of latticeconstants of the crystal in the three directions. • transverse-truncation determines a limit to transversely truncate the bunches. This factor brings the possibilityto control particle initialization and prevents them from escaping out of the computational domain. The bunch initializertruncates the bunch at the given distance from the bunch center. • longitudinal-truncation determines a limit to longitudinally truncate the bunches. This factor brings the pos-sibility to control particle initialization and prevents them from escaping out of the computational domain. The bunchinitializer truncates the bunch at the given distance from the bunch center. • bunching-factor is a value larger than zero and less than one, which determines the bunching factor, i.e. < e jk u z > ,of the initialized bunch.2. bunch-sampling: This group defines the required parameters for saving the bunch properties with time. The bunchmean position, mean momentum, position spread, and momentum spread along the three Cartesian coordinates are savedrespectively with time. There are different parameters required for this definition which include: • sample is a boolean value determining if the bunch sampling should be activated. • base-name is the file name (no suffix) with the required address of the file to save the output data. • directory is the address where the above file should be saved. The file name with the address can also be given in thebase-name section. The software eventually considers the combination of directory and base-name as the final completefile name. • rhythm is the rhythm of bunch sampling, i.e. the time interval between two consecutive visualization times.3. bunch-visualization: This group defines the required parameters for visualizing the charge distribution in the wholecomputational domain. The output will be a set of .vtu files at each time for each processor which are connected with a setof .pvtu files. They can be very nicely visualized using the open source ParaView package. There are different parametersrequired for this definition which include: • sample is a boolean value determining if the charge visualization should be activated. • base-name is the file name (no suffix) with the required address of the file to save the output data. • directory is the address where the above file should be saved. The file name with the address can also be given in thebase-name section. The software eventually considers the combination of directory and base-name as the final completefile name. • rhythm is the rhythm of charge illustration, i.e. the time interval between two consecutive visualization times.4. bunch-profile: This group defines the required parameters for saving a histogram of the charges. It means that at aspecific time instant the charge values, positions and momenta of all the particles (or macro-particles) will be written andsaved in a file. The parameters entered by the user for saving the histogram include: • sample is a boolean value determining if the writing of the histogram during the PIC simulations should be activated. • base-name is the file name (no suffix) with the required address of the file to save the output data. • directory is the address where the above file should be saved. The file name with the address can also be given in thebase-name section. The software eventually considers the combination of directory and base-name as the final completefile name. • time is the time instant for saving the histogram. If this needs to be done in several time instants, simply this line shouldbe repeated with different time values. HAPTER 3. USER INTERFACE 25 • rhythm is the rhythm of writing the bunch profile, i.e. the time interval between two consecutive profiling times. Ifthis value is nonzero, the sequence of times will be considered in addition to the specific time points given by the timevariable.The format of the BUNCH group is:
BUNCH{ bunch-initialization{ type = < manual | ellipsoid | 3D-crystal | file >distribution = < uniform | gaussian >charge = < real value >number-of-particles = < integer value >gamma = < real value >beta = < real value >direction = < ( real value, real value, real value ) >position = < ( real value, real value, real value ) >sigma-position = < ( real value, real value, real value ) >sigma-momentum = < ( real value, real value, real value ) >transverse-truncation = < real value >longitudinal-truncation = < real value >bunching-factor = < real value between zero and one >}bunch-sampling{ sample = < true | false >directory = < address according to UNIX convention >base-name = < name of the file >rhythm = < real value >}bunch-visualization{ sample = < true | false >directory = < address according to UNIX convention >base-name = < name of the file >rhythm = < real value >}bunch-profile{ sample = < true | false >directory = < address according to UNIX convention >base-name = < name of the file >time = < real value >rhythm = < real value >}}
An example of the bunch category definition looks as the following:
BUNCH{ bunch-initialization{ type = ellipsoiddistribution = uniform
HAPTER 3. USER INTERFACE 26 charge = 1.846e8number-of-particles = 131072gamma = 100.41direction = ( 0.0, 0.0, 1.0)position = ( 0.0, 0.0, 0.0)sigma-position = ( 260.0, 260.0, 50.25)sigma-momentum = ( 1.0e-8, 1.0e-8, 100.41e-4)transverse-truncation = 1040.0longitudinal-truncation = 90.0bunching-factor = 0.01}bunch-sampling{ sample = falsedirectory = ./base-name = bunch-sampling/bunchrhythm = 3.2}bunch-visualization{ sample = truedirectory = ./base-name = bunch-visualization/bunchrhythm = 32}bunch-profile{ sample = falsedirectory = ./base-name = bunch-profile/bunchtime = 5000time = 10000time = 15000time = 20000time = 25000time = 30000}}
FIELD
In section
FIELD , the required data for the input and output framework of the field in the FDTD algorithm is produced. Thissection consists of four groups: (1) field-initialization , (2) field-sampling , (3) field-visualization , and(4) field-profile . As apparent from the name, the first group determines the set of parameters to initialize the field and theother three groups are dedicated to reporting the field propagation in different formats. In what follows, the parameters in each groupare introduced:1. field-initialization:
This group mainly determines the parameters whose values are needed for initializing a fieldexcitation entering the computational domain. The excitation may have different types. This group is where a seed can beadded to the simulations to simulate a seeded-FEL problem. The set of values accepted in this group include: • type is the type of the excitation. The accepted excitation types in MITHRA include plane wave and Gaussian beam. • position is the reference position of the excitation. It is the reference position of the plane wave propagation in theplane wave excitation and the focusing point in the Gaussian beam excitation. HAPTER 3. USER INTERFACE 27 • direction is the propagation direction of the excitation in the plane wave and Gaussian beam types. • polarization is the polarization of the incoming excitation and is used by both plane wave and Gaussian beam types. • rayleigh-radius-parallel is the Rayleigh radius of the Gaussian beam in the direction parallel to the polariza-tion and is only meaningful for the Gaussian beam type. • rayleigh-radius-perpendicular is the Rayleigh radius of the Gaussian beam in the direction perpendicularto the polarization and is only meaningful for the Gaussian beam type. • signal-type determines the time signature of the signal exciting the fields according to the particular type. Theaccepted signal types in MITHRA include modulated Neumann, modulated Gaussian, modulated secant hyperbolic andthe sinusoidal pulse. The equation representing the time domain variation of each pulse is as follows:modulated Neumann: f ( t ) = − A πf ( t − t ) + φ CEP ) t − t τ e − t − t ) /τ modulated Gaussian: f ( t ) = A cos(2 πf ( t − t ) + φ CEP ) e − t − t ) /τ modulated hyperbolic secant: f ( t ) = A cos(2 πf ( t − t ) + φ CEP ) 1cosh(( t − t ) /τ ) sinusoidal pulse: f ( t ) = A cos(2 πf ( t − t ) + φ CEP ) e − t − t ) /τ t ≤ t A cos(2 πf ( t − t ) + φ CEP ) t > t (3.1) • strength-parameter is the normalized amplitude a = eA /m e c of the beam. • offset is the distance offset of the signal ct . • variance is the variance of the signal in length units cτ . • wavelength is the modulation wavelength λ of the modulated signal. • CEP is the carrier envelope phase φ CEP of the modulated signal.2. field-sampling:
This group defines the required parameters for saving the field value at specific points with time. Thereare different parameters required for this definition which include: • sample is a boolean value determining if the field sampling should be activated. • type determines if the field should be sampled at the given points (at-point) or the field should be sampled at the pointsover a line (over-line). • field determines which electromagnetic field is to be sampled. The available options are the electric field, magneticfield, magnetic vector potential, scalar electric potential, charge and current. This item can be repeated to assign severalfields for the sampling. In the text file, the fields appear with the same order. • base-name is the file name (no suffix) with the required address of the file to save the output data. • directory is the address where the above file should be saved. The file name with the address can also be given in thebase-name section. The software eventually considers the combination of directory and base-name as the final completefile name. • rhythm is the rhythm of field sampling, i.e. the time interval between two consecutive visualization times. • position is the coordinate of the points where the fields should be sampled. By repeating this line any numberof points can be added to the set of sampling locations. This option is merely used when the sampling type is set to”at-point” option. • line-begin defines the position of the line begin over which the fields should be sampled and is used when thesampling type is set to ”over-line” option. • line-end defines the position of the line end over which the fields should be sampled and is used when the samplingtype is set to ”over-line” option. • resolution defines the resolution of the line discretization for sampling. In other words, the value is the distancebetween two adjacent sampling points. This value is used when the sampling type is set to ”over-line” option. HAPTER 3. USER INTERFACE 28 field-visualization: This group defines the required parameters for visualizing the fields in the whole computationaldomain. The output will be a set of .vtu files at each time for each processor which are connected with a set of .pvtu files.They can be very nicely visualized using the open source ParaView package. There are different parameters required for thisdefinition which include: • sample is a boolean value determining if the field visualization should be activated. • base-name is the file name (no suffix) with the required address of the file to save the output data. • directory is the address where the above file should be saved. The file name with the address can also be given in thebase-name section. The software eventually considers the combination of directory and base-name as the final completefile name. • rhythm is the rhythm of field visualization, i.e. the time interval between two consecutive visualization times. • field determines which electromagnetic field is to be saved. The available options are the electric field, magnetic field,magnetic vector potential, scalar electric potential, charge and current. This item can be repeated to assign several fieldsfor the sampling. In the text file, the fields appear with the same order.4. field-profile: This group defines the required parameters for saving a histogram of the field over the whole computa-tional domain. It means that at a specific time instant the field values and the corresponding positions at all the grid points willbe written and saved in a text file. The parameters entered by the user for saving the histogram include: • sample is a boolean value determining if the writing of the histogram during the FDTD simulations should be activated. • base-name is the file name (no suffix) with the required address of the file to save the output data. • directory is the address where the above file should be saved. The file name with the address can also be given in thebase-name section. The software eventually considers the combination of directory and base-name as the final completefile name. • time is the time instant for saving the histogram. If this needs to be done in several time instants, simply this line shouldbe repeated with different time values. • rhythm is the rhythm of field profiling, i.e. the time interval between two consecutive visualization times. Bothrhythmic profiling and saving the fields at specific times can be given to the software. • field determines which electromagnetic field is to be profiled. The available options are the electric field, magneticfield, magnetic vector potential, scalar electric potential, charge and current. This item can be repeated to assign severalfields for the sampling. In the text file, the fields appear with the same order.The format of the FIELD group is:
FIELD{ field-initialization{ type = < hertzian-dipole | plane-wave | gaussian-beam >position = < ( real value , real value , real value ) >direction = < ( real value , real value , real value ) >polarization = < ( real value , real value , real value ) >rayleigh-radius-parallel = < real value >rayleigh-radius-perpendicular = < real value >signal-type = < neumann | gaussian | secant-hyperbolic | flat-top >strength-parameter = < real value >offset = < real value >variance = < real value >wavelength = < real value >CEP = < real value >}field-sampling{ sample = < true | false >
HAPTER 3. USER INTERFACE 29 type = < over-line | at-point >field = < Ex | Ey | Ez | Bx | By | Bz | Ax | Ay | Az | Jx | Jy| Jz | F | Q >directory = < address according to UNIX convention >base-name = < name of the file >rhythm = < real value >position = < ( real value , real value , real value ) >line-begin = < ( real value , real value , real value ) >line-end = < ( real value , real value , real value ) >resolution = < real value >}field-visualization{ sample = < true | false >field = < Ex | Ey | Ez | Bx | By | Bz | Ax | Ay | Az | Jx | Jy| Jz | F | Q >directory = < address according to UNIX convention >base-name = < name of the file >rhythm = < real value >}field-profile{ sample = < true | false >field = < Ex | Ey | Ez | Bx | By | Bz | Ax | Ay | Az | Jx | Jy| Jz | F | Q >directory = < address according to UNIX convention >base-name = < name of the file >rhythm = < real value >time = < real value >}}
An example of the field category definition looks as the following:
FIELD{ field-initialization{ type = gaussian-beamposition = ( 0.0, 0.0, -2500.0)direction = ( 0.0, 0.0, 1.0)polarization = ( 0.0, 1.0, 0.0)rayleigh-radius-parallel = 0.5rayleigh-radius-perpendicular = 0.5strength-parameter = 0.0signal-type = gaussianoffset = 0.00variance = 1.00wavelength = 0.0CEP = 0.0}field-sampling{ sample = truetype = at-pointfield = Ex
HAPTER 3. USER INTERFACE 30 field = Eyfield = Ezdirectory = ./base-name = field-sampling/fieldrhythm = 3.2position = (0.0, 0.0, 110.0)}field-visualization{ sample = truefield = Exfield = Eyfield = Ezfield = Qdirectory = ./base-name = field-visualization/fieldrhythm = 32}field-profile{ sample = falsefield = Exfield = Eyfield = Ezdirectory = ./base-name = field-profile/fieldrhythm = 100}}
UNDULATOR
In section
UNDULATOR , the properties of the undulator considered in the FEL problem are introduced. The parameters, throughwhich the complete data for establishing undulator fields are obtained, include: • undulator-type is the type of the undulator considered for the FEL interaction. It can be either an optical or staticundulator. • undulator-parameter is the undulator parameter of the undulator, i.e. the so-called K parameter. • period is the period of the undulator in the given length-scale determined in the mesh class. • length is the total length of the undulator. • polarization-angle is the angle between the magnetic field polarization and the x -axis in degrees. • beam-type is the type of the pulse for an optical undulator. The accepted undulator beam types in MITHRA include planewave and Gaussian beam. • position is the reference position of the undulator. It is the reference position of the plane wave propagation in the planewave undulator and the focusing point in the Gaussian beam undulator. • direction is the propagation direction of the optical undulator in the plane wave and Gaussian beam types. • polarization is the polarization of the undulator and is used by both plane wave and Gaussian beam types. • rayleigh-radius-parallel is the Rayleigh radius of the Gaussian beam in the direction parallel to the polarizationand is only meaningful for the Gaussian beam type. HAPTER 3. USER INTERFACE 31 • rayleigh-radius-perpendicular is the Rayleigh radius of the Gaussian beam in the direction perpendicular to thepolarization and is only meaningful for the Gaussian beam type. • signal-type determines the time signature of the undulator. The accepted signal types in MITHRA are listed in the fieldsection. Here, the same set of signal can be given as an undulator envelope. • strength-parameter is the normalized amplitude a = eA /m e c of the undulator, which is equivalent to the undulator-parameter in the static case. • offset is the distance offset of the signal ct . • variance is the time variance of the signal in length unit cτ . • wavelength is the modulation wavelength λ of the modulated signal. • CEP is the carrier envelope phase φ CEP of the modulated signal.The format of the
UNDULATOR group is:
UNDULATOR{ undulator-type = < static | optical >undulator-parameter = < real value >period = < real value >length = < real value >polarization-angle = < real value >beam-type = < hertzian-dipole | plane-wave | gaussian-beam >position = < ( real value , real value , real value ) >direction = < ( real value , real value , real value ) >polarization = < ( real value , real value , real value ) >rayleigh-radius-parallel = < real value >rayleigh-radius-perpendicular = < real value >signal-type = < neumann | gaussian | secant-hyperbolic | flat-top >strength-parameter = < real value >offset = < real value >variance = < real value >wavelength = < real value >CEP = < real value >}
As explained before, MITHRA always initializes the bunch outside the undulator. It may be already noticed that there exists nooption to determine the beginning of the static undulator with respect to the bunch. This is automatically set by the solver, to avoidparticle initialization inside the undulator. For the optical undulator type, the user should control this effect through the parameteroffset. An example of the undulator category definition looks as the following:
UNDULATOR{ undulator-type = staticundulator-parameter = 1.417period = 3.0e4length = 300polarization-angle = 0.0}
An instance of optical undulator definition reads as follows:
UNDULATOR{ undulator-type = optical
HAPTER 3. USER INTERFACE 32 beam-type = plane-waveposition = ( 0.0, 0.0, 0.0 )direction = ( 0.0, 0.0,-1.0 )polarization = ( 0.0, 1.0, 0.0 )strength-parameter = 0.5signal-type = flat-topwavelength = 1.0e3variance = 1200.0e3offset = 600118.0CEP = 0.0}
FEL-OUTPUT
The typical parameters for a free electron laser instrument are calculated from the radiated fields using the parameter definitionsat this section. Currently, there is only one group implemented in MITHRA, which calculates the total radiated power in the + z direction and saves the power versus the undulator position to a text file. In what follows, the parameters in each group are introduced:1. radiation-power: This group mainly determines the parameters whose values are needed for calculating the radiationpower from the field distribution. The output is a .txt file with the first column the time point and the second column theradiated power at the given point. For multiple points the column is repeated. If multiple wavelengths are given there willbe several rows with similar time value listing the radiated power with different wavelengths. This group can be repeated toobtain various files for different output definitions.The set of values accepted in this group include: • sample is a boolean parameter which activates the computation of the total radiated power at each instant. • base-name is the file name (no suffix) with the required address of the file to save the output data. • directory is the address where the above file should be saved. The file name with the address can also be given in thebase-name section. The software eventually considers the combination of directory and base-name as the final completefile name. • type determines if the radiated power should be sampled at given distances from the bunch (at-point) or should besampled at the points over a line (over-line). • distance-from-bunch gives the total distance from the bunch where a sampling plate to capture the whole radiatedpower will be placed. One can enter several sampling positions by repeating this line. This option is only considered ifthe type variable is set to ”at-point”. • line-begin defines the distance from the bunch center for the line begin over which the fields should be sampled andis used when the sampling type is set to ”over-line” option. In case of multiple bunches, the center of the first bunch isconsidered as the reference position. • line-end defines the distance from the bunch center for the line end over which the fields should be sampled and isused when the sampling type is set to ”over-line” option. In case of multiple bunches, the center of the first bunch isconsidered as the reference position. • resolution defines the resolution of the line discretization for sampling. In other words, the value is the distancebetween two adjacent sampling points. This value is used when the type variable is set to ”over-line” option. • normalized-frequency is the central wavelength normalized to the radiation wavelength of the radiation spectrum. • minimum-normalized-frequency is the minimum wavelength normalized to the radiation wavelength. Thisparameters and the next two parameters are used to sweep over the normalized wavelength and save the radiation powerspectrum. • maximum-normalized-frequency is the maximum wavelength normalized to the radiation wavelength. • normalized-frequency-resolution is the sweep resolution for the wavelength normalized to the radiationwavelength. HAPTER 3. USER INTERFACE 33
The format of the
FEL-OUTPUT group is:
FEL-OUTPUT{ radiation-power{ sample = < false | true >type = < at-point | over-line >directory = < address according to UNIX convention >base-name = < name of the file >distance-from-bunch = < real value >line-begin = < real value >line-end = < real value >resolution = < real value >normalized-frequency = < real value >minimum-normalized-frequency = < real value >maximum-normalized-frequency = < real value >normalized-frequency-resolution = < real value >}}
An example of the FEL output category definition looks as the following:
FEL-OUTPUT{ radiation-power{ sample = truetype = at-pointdirectory = ./base-name = power-sampling/powerdistance-from-bunch = 110.0normalized-frequency = 1.00}} hapter 4
Examples
The goal in this chapter is to present several examples for the MITHRA users to more easily get familiar with the interface of thesoftware. In addition, through the presented examples the pros and cons of using the developed FDTD/PIC algorithm are moreaccurately evaluated. For example, the computation time, numerical stability and numerical convergence and more importantly thereliability of the results are studied based on some standard examples. The software developers aim to update this chapter with themost recent examples where MITHRA is used for the FEL simulation.
Table 4.1: Parameters of the Infrared FEL configuration considered as the first example.FEL parameter ValueCurrent profile UniformBunch size (260 × × µ mBunch charge 29.5 pCBunch energy 51.4 MeVBunch current 88.5 ALongitudinal momentum spread 0.01%Transverse momentum spread 0.0Undulator period 3.0 cmMagnetic field 0.5 TUndulator parameter 1.4Undulator length 5 mRadiation wavelength 3 µ mElectron density . × / cm Gain length (1D) 22.4 cmFEL parameter 0.006Cooperation length . µ mInitial bunching factor . As the first example, we consider an infrared FEL with the parameters tabulated in table 4.1.1, which is inspired by the numericalanalysis presented in [9]. The bunch distribution is assumed to be uniform in order to compare the results with one-dimensionalFEL theory. For the same purpose, the transverse energy spread is considered to be zero and a minimal longitudinal energy spreadis assumed.To simulate the considered FEL configuration, the following job file is written and given to the software to analyze the interaction.As observed in the mesh definition, the transverse size of the computational domain is almost 10 times larger than the bunchtransverse size. In the contrary, the longitudinal size of the mesh is only three times larger than the bunch length. This needs to34
HAPTER 4. EXAMPLES 35 be considered due to the failure of absorbing boundary conditions for the oblique incidence of the field. Furthermore, the bunchand undulator both have some tapering sections to avoid abrupt transitions which produce coherent scattering emission (CSE). It isrecommended that the undulator begin is initialized at least ten radiation wavelengths apart from the bunch head to reduce the CSE.This also introduces corresponding limitations on the mesh size, meaning that the minimum distance from the bunch tail and themesh boundary should be at least ten radiation wavelengths. In the illustrated job file, some of the output formats are turned offwhich can always be activated to obtain the required data.
MESH{ length-scale = MICROMETERtime-scale = PICOSECONDmesh-lengths = ( 3200, 3200.0, 280.0)mesh-resolution = ( 50.0, 50.0, 0.1)mesh-center = ( 0.0, 0.0, 0.0)total-time = 30000bunch-time-step = 1.6bunch-time-start = 0.0number-of-threads = 2mesh-truncation-order = 2space-charge = false}BUNCH{ bunch-initialization{ type = ellipsoiddistribution = uniformcharge = 1.846e8number-of-particles = 131072gamma = 100.41direction = ( 0.0, 0.0, 1.0)position = ( 0.0, 0.0, 0.0)sigma-position = ( 260.0, 260.0, 50.25)sigma-momentum = ( 1.0e-8, 1.0e-8, 100.41e-4)transverse-truncation = 1040.0longitudinal-truncation = 90.0bunching-factor = 0.01}bunch-sampling{ sample = falsedirectory = ./base-name = bunch-sampling/bunchrhythm = 0.01}bunch-visualization{ sample = truedirectory = ./base-name = bunch-visualization/bunchrhythm = 32}bunch-profile{ sample = false
HAPTER 4. EXAMPLES 36 directory = ./base-name = bunch-profile/bunchtime = 15.0}}FIELD{ field-initialization{ type = gaussian-beamposition = ( 0.0, 0.0, -2500.0)direction = ( 0.0, 0.0, 1.0)polarization = ( 0.0, 1.0, 0.0)rayleigh-radius-parallel = 0.5rayleigh-radius-perpendicular = 0.5strength-parameter = 0.0signal-type = gaussianoffset = 0.00e-9variance = 1.00e-100wavelength = 0.0CEP = 0.0}field-sampling{ sample = truetype = at-pointfield = Exfield = Eyfield = Ezdirectory = ./base-name = field-sampling/fieldrhythm = 3.2position = (0.0, 0.0, 110.0)}field-visualization{ sample = truefield = Exfield = Eyfield = Ezfield = Qdirectory = ./base-name = field-visualization/fieldrhythm = 32}field-profile{ sample = falsefield = Exfield = Eyfield = Ezdirectory = ./base-name = field-profile/fieldrhythm = 100}
HAPTER 4. EXAMPLES 37 }UNDULATOR{ undulator-parameter = 1.417period = 3.0e4length = 300polarization-angle = 0.0}FEL-OUTPUT{ radiation-power{ sample = falsetype = at-pointdirectory = ./base-name = power-sampling/powerdistance-from-bunch = 110.0normalized-frequency = 1.00}}
In the beginning, we neglect the space-charge effect only to achieve a good assessment of MITHRA simulation results. The investi-gation of space-charge effect will be performed in the second step. Fig. 4.1a shows the transverse electric field sampled at 110 µ min front of the bunch center. The logarithmic plot of the radiated power for different propagation lengths ( z ) is also depicted inFig. 4.1b. According to the 1D FEL theory the gain length of the considered SASE FEL configuration is L G = 22 . . The gainlength calculated from the slope of the power curve is L G = 22 cm . There exists also a good agreement in the computed satura-tion power. The beam energy according to the data in table 4.1.1 is 1.52 mJ which for the bunch length of 100 µ m corresponds to P beam = 4 .
55 GW beam power. The estimated saturation power according to the 1D theory is equal to P sat = ρP beam = 2 . .The saturation power computed by MITHRA is . .We have also performed a comparison study between the obtained results from MITHRA and the code Genesis 1.3, which ispresented in Fig. 4.1b. As observed, both codes produce similar results in the initial state and the gain regime. Nonetheless, thereexists a considerable discrepancy between the calculated radiated power in the saturation regime. The illustrated results in Fig. 4.1bshow that the steady state and time domain analyses using Genesis do not produce similar results. This shows that the bunch is notlong enough to justify the steady state approximation, and dictates a time domain analysis for accurate simulation. However, theresults obtained by MITHRA at saturation do not match with the Genesis results even in the time domain.The origin of such a discrepancy is described as follows: As explained in chapter 1, Genesis 1.3 and all the existing softwaresfor FEL simulation neglect the backward radiation of the electrons. Such an approximation is motivated by the inherent interestin forward radiation in the FEL process. However, the backward radiation although is seldom used due to its long wavelength, itinfluences the motion of electrons, the charge distribution and in turn the FEL output. The influence of low-frequency backwardradiation on the performance of free electron lasers has been already studied in a 1D regime [41]. The effect becomes strongerin the saturation regime, when the electrons are microbunched and the FEL forward radiation is a strong function of the particlesdistribution. To demonstrate this effect, we used a mathematical trick in MITHRA through the domain decomposition algorithm tosuppress the propagation of backward radiation. The results of such an analysis is also shown in Fig. 4.1b, which shows a relativelybetter agreement with time domain simulation results returned by Genesis 1.3. The still existing discrepancy is attributed to thedifferent formulations of FDTD and TDA algorithms as well as the introduced tapers in bunch current and undulator fields.As a 3D electromagnetic simulation, it is always beneficial to investigate the electromagnetic field profile in the computationaldomain. Using the field visualization capability in MITHRA, snapshots of the field profile at different instants and from variousview points are provided. In Fig. 4.2, snapshots of the radiated field profile at different time instants are illustrated. The emergenceof lasing radiation at the end of the undulator motion is clearly observed in the field profile.Furthermore, snapshots of the bunch profile are also presented beside the field profile. The main FEL principle which is thelasing due to micro-bunching of the electron bunch is observed from the field and bunch profiles. HAPTER 4. EXAMPLES 38 z (m) -50050 E z ( M V / m ) z rm) l og P [ W ] Genesis 1.3 time domainGenesis 1.3 steady stateMITHRAMITHRA BW suppressed ( a ) ( b ) Figure 4.1: (a) The transverse field E y at 110 µ m distance from the bunch center and (b) the total radiated power calculated at110 µ m distance from the bunch center in terms of the traveled undulator length. (a) (b) (c)
20 mm 3.2 mm 10 mm
Figure 4.2: Snapshots of the radiated field profile taken at (a) x = 0 and (b) z = 60 µm plane and (c) the bunch profile viewed fromthe x axis. HAPTER 4. EXAMPLES 39 l x = l y ( µ m) e rr o r( % ) ∆ x=1/ ∆ y (1/ µ m) e rr o r( % ) n × e rr o r( % )
10 20 30 40 ∆ z (1/ µ m) e rr o r( % ) ∆ t b (1/fs) e rr o r( % ) (a)(d) (b) (c)(e) (f)
230 240 250 260 270 280 l z ( µ m) e rr o r( % ) (c) Figure 4.3: Convergence study for the different involved parameters in the considered FEL simulation: (a) n , (b) ∆ t b , (c) l z , (d) l x = l y , (e) ∆ z and (f) ∆ x = ∆ y The convergence rate of the results is the main factor used to assess a numerical algorithm. In our FEL analysis, there are severalparameters introduced by the numerical method which may affect the final result. These parameters include (1) number of macro-particles ( n ), (2) time step for updating equation of motion ( ∆ t b ), (3) longitudinal mesh size ( l z ), (4) transverse mesh size ( l x = l y ),(5) longitudinal discretization ( ∆ z ) and (6) transverse discretization ( ∆ x = ∆ y ). Studying the convergence of the results is crucial toacquire an estimate for the uncertainty in the reported values. Here, this task is accomplished by sweeping over the above parametersand plotting the error function defined as the following: error = ´ z f z i | P ( z ) − P ( z ) | dz ´ z f z i P ( z )dz , (4.1)where z i and z f are the beginning and end of the undulator, respectively and P is the reference simulation result which is chosen asthe results with the highest resolution.In Fig. 4.3 the convergence study is shown for the aforementioned parameters. Generally, accuracy of less than 3% is achievedby using the initially suggested values. A promising benefit offered by MITHRA is the assessment of various approximations used in the previously developed FEL codes.As an example, the algorithm used in the TDA method to evaluate the space-charge effect can be examined and verified using thiscode. The TDA method implemented in Genesis 1.3 software considers a periodic variation of space-charge force throughout theelectron bunch [39, 42]. However, a simple investigation of bunch profiles shown in Fig. 4.2c shows that a periodic assumptionfor the electron distribution may be a crude approximation. In addition, this assumption is favored by the FEL gain process andpotentially decreases any detrimental influence of the space-charge fields on the FEL radiation. In Fig. 4.4a, we are comparing thesolution of the FEL problem using MITHRA and Genesis 1.3 with and without considering the space-charge effect. As observed inthe results, the effect of space-charge on the radiation gain predicted by MITHRA is much stronger than the same effect predictedby Genesis. This is attributed to the assumption of periodic variations in the space-charge force made in TDA algorithm. If such
HAPTER 4. EXAMPLES 40 l og P [ W ] Genesis 1.3 no SCGenesis 1.3 with SCMITHRA no SCMITHRA with SC l og P [ W ] Genesis 100 µ mGenesis 200 µ mGenesis 400 µ mMITHRA 100 µ mMITHRA 200 µ mMITHRA 400 µ m ( a ) ( b ) Figure 4.4: The total radiated power calculated at 110 µ m distance from the bunch center in terms of the traveled undulator length(a) with and without space-charge consideration and (b) various lengths of the bunch with space-charge assumption.a hypothesis is correct, the observed discrepancy should reduce once the radiation from a longer bunch is simulated, because theaccuracy of periodicity assumption increases for longer bunches. Indeed, this is observed after repeating the simulation for longerelectron bunches with similar charge and current densities. The results of such a study is illustrated in Fig. 4.4b. As the second example, we consider a seeded FEL in the UV regime to verify the implemented features for simulating a seededFEL. The parameters of the considered case are taken from [12], which are tabulated in table 4.2.1. The bunch distribution is againassumed to be uniform with a long current profile ( ∼ MESH{ length-scale = MICROMETERtime-scale = PICOSECONDmesh-lengths = ( 1600.0, 1600.0, 165.0)mesh-resolution = ( 10.0, 10.0, 0.01 )mesh-center = ( 0.0, 0.0, 0.0 )total-time = 50000bunch-time-step = 0.8bunch-time-start = 0.0number-of-threads = 2mesh-truncation-order = 2space-charge = true}
HAPTER 4. EXAMPLES 41
Table 4.2: Parameters of the UV seeded FEL configuration considered as the second example.FEL parameter ValueCurrent profile UniformBunch size (95 . × . × µ mBunch charge 54.9 pCBunch energy 200 MeVBunch current 110 ALongitudinal momentum spread 0.01%Transverse momentum spread 0.0026%Undulator period 2.8 cmMagnetic field 0.7 TUndulator parameter 1.95Undulator length 15 mRadiation wavelength 0.265 µ mElectron density . × / cm Gain length (1D) 8.9 cmFEL parameter 0.015Cooperation length . µ mInitial bunching factor . Seed type Gaussian beamSeed focal point 70 cmSeed beam radius 183.74 µ mSeed pulse length infiniteSeed power 10 kW BUNCH{ bunch-initialization{ type = ellipsoiddistribution = uniformcharge = 3.4332e8number-of-particles = 1048576gamma = 391.36direction = ( 0.0, 0.0, 1.0 )position = ( 0.0, 0.0, 0.0 )sigma-position = ( 95.3, 95.3, 75.0)sigma-momentum = ( 0.0105, 0.0105, 391.36e-4)transverse-truncation = 400.0longitudinal-truncation = 78.0bunching-factor = 0.0}bunch-sampling{ sample = falsedirectory = ./base-name = bunch-sampling/bunchrhythm = 0.01}bunch-visualization{ sample = false
HAPTER 4. EXAMPLES 42 directory = ./base-name = bunch-visualization/bunchrhythm = 10}bunch-profile{ sample = truedirectory = ./base-name = bunch-profile/bunchtime = 40000}}FIELD{ field-initialization{ type = gaussian-beamposition = ( 0.0, 0.0, 700080)direction = ( 0.0, 0.0, 1.0)polarization = ( 0.0, 1.0, 0.0)rayleigh-radius-parallel = 183.74rayleigh-radius-perpendicular = 183.74strength-parameter = 9.857e-7signal-type = gaussianoffset = 700005.0variance = 1e12wavelength = 0.265187CEP = 0.0}field-sampling{ sample = falsetype = at-pointfield = Exdirectory = ./base-name = field-sampling/fieldrhythm = 0.5position = (0.0, 0.0, 55.0)}field-visualization{ sample = falsefield = Aydirectory = ./base-name = field-visualization/fieldrhythm = 100}field-profile{ sample = falsefield = Exdirectory = ./base-name = field-profile/fieldrhythm = 100
HAPTER 4. EXAMPLES 43 z (m) l og P [ W ] Genesis no SCGenesis with SCMITHRA no SCMITHRA with SC ( a ) ( b ) Figure 4.5: (a) The total radiated power measured at 80 µ m distance from the bunch center in terms of the traveled undulator lengthand (b) the bunch profile at 12 m from the undulator begin. }}UNDULATOR{ undulator-parameter = 1.95period = 2.8e4length = 535polarization-angle = 0.0}FEL-OUTPUT{ radiation-power{ sample = truetype = at-pointdirectory = ./base-name = power-sampling/powerdistance-from-bunch = 80.0normalized-frequency = 1.00}} Fig. 4.5a shows the radiated power in terms of travelled undulator distance computed using MITHRA and Genesis. As observedagain in this example, the results agree very well in the seeded and gain regime, with notable discrepancies in the saturation regime.In Fig. 4.5b, the bunch profile after 12 m propagation in the undulator is also depicted. The micro-bunching of the large bunch is onlyvisible once a zoom into a part of the bunch is considered. The investigation of the results with and without considering space-chargeeffect shows that in the seeded and gain intervals, space charge plays a negligible role. However, in the saturation regime the effectof space-charge predicted by MITHRA is stronger than the effect predicted by Genesis.
HAPTER 4. EXAMPLES 44 z (m) l og P [ W ] Static undulatorOptical undulator
Figure 4.6: The total radiated power calculated at 110 µ m distance from the bunch center in terms of the traveled undulator lengthcompared for two cases of an optical and static undulator. As explained in the introduction of this manual, one of the milestones considered for the development of MITHRA is full-wavesimulation of inverse Compton scattering (ICS) or the so-called optical undulator. The possibility of lasing or the so-called micro-bunching in an electron beam due to an interaction with a counter-propagating laser beam has been under debate for several years. Afull-wave analysis of such an interaction definitely gives valuable physical insight to this process. Note that the classical treatment ofthis interaction within MITHRA does not allow for any consideration of quantum mechanical effects. It is known that the radiationof photons results in a backward force on electrons which leads to a change in their momenta. This effect is considered in thesimulations presented by MITHRA within a classical framework. Nonetheless, the fact that this momentum change is a quantizedvalue and equal to n (cid:126) ω can not be considered by MITHRA.Before embarking on the analysis and interpretation of results for a typical ICS experiment, a benchmark to validate the analysisof optical undulators using FDTD/PIC is presented. It is known that electron trajectories in a static undulator with undulator pa-rameter K and periodicity λ u are similar to the trajectories in an electromagnetic undulator setup with normalized vector potential a = K and wavelength λ l = 2 λ u [43]. We take the first SASE FEL example in table 4.1.1 into account and analyze the sameconfiguration but with an equivalent optical undulator. For this purpose, the undulator definition of example 1 is entered accordingto the following: UNDULATOR{ undulator-type = opticalbeam-type = plane-waveposition = ( 0.0, 0.0, 0.0 )direction = ( 0.0, 0.0,-1.0 )polarization = ( 0.0, 1.0, 0.0 )strength-parameter = 1.417signal-type = flat-topwavelength = 6.0e4variance = 1800.0e4offset = 9000150CEP = 0.0}
Fig. 4.6 illustrates a comparison between the radiated infrared light for the static and optical undulator cases. The very closeagreement between the two results validates the implementation of optical undulators in MITHRA.The parameters of FEL interaction in an optical undulator, considered as the third example, are tabulated in table 4.3. Sincewe observe drastic deviation from the predictions of one-dimensional FEL theory in our simulations, we have not listed the FEL
HAPTER 4. EXAMPLES 45 parameters calculated using the 1D theory. We believe the discrepancies are originated from the small number of electrons in each3D wave bucket, i.e. only 2 electrons. This strongly intensifies the 3D effects, dramatically reduces the transverse coherence ofthe radiation, and indeed makes analysis using 1D FEL theory completely invalid. Detailed verification of the above claims will bepresented in a publication. Here, the optical undulator is merely treated as an example analyzed with FDTD/PIC.Table 4.3: Parameters of the FEL configuration with optical undulator considered as the third example.FEL parameter ValueCurrent profile UniformBunch size (60 × × nmBunch charge 0.45 fCBunch energy 15 MeVBunch current 0.93 ALongitudinal momentum spread 0.1%Transverse momentum spread 0.1%Laser wavelength 1 µ mLaser strength parameter 1.0Pulse duration 4 psLaser pulse type flat-topRadiation wavelength 0.41 nmElectron density . × / cm Initial bunching factor . To simulate the considered FEL configuration, the following job file is written and given to the software to analyze the interaction.
MESH{ length-scale = NANOMETERtime-scale = ATTOSECONDmesh-lengths = ( 2000.0, 2000.0, 165.0 )mesh-resolution = ( 5.0, 5.0, 0.05)mesh-center = ( 0.0, 0.0, 0.0 )total-time = 2000000bunch-time-step = 20.0bunch-time-start = 0.0number-of-threads = 4mesh-truncation-order = 2space-charge = true}BUNCH{ bunch-initialization{ type = ellipsoiddistribution = uniformcharge = 2800number-of-particles = 2800gamma = 30.0direction = ( 0.0, 0.0, 1.0)position = ( 0.0, 0.0, 0.0)sigma-position = ( 60.0, 60.0, 72.0)sigma-momentum = ( 0.03, 0.03, 0.03)transverse-truncation = 240.0longitudinal-truncation = 77.0bunching-factor = 0.0}
HAPTER 4. EXAMPLES 46 bunch-sampling{ sample = falsedirectory = ./base-name = bunch-sampling/bunchrhythm = 1000}bunch-visualization{ sample = falsedirectory = ./base-name = bunch-visualization/bunchrhythm = 1600}bunch-profile{ sample = falsedirectory = ./base-name = bunch-profile/bunchrhythm = 1600}}FIELD{ field-initialization{ type = gaussian-beamposition = ( 0.0, 0.0, -2500.0)direction = ( 0.0, 0.0, 1.0)polarization = ( 0.0, 1.0, 0.0)rayleigh-radius-parallel = 0.5rayleigh-radius-perpendicular = 0.5strength-parameter = 0.0signal-type = gaussianoffset = 0.00e-9variance = 1.00wavelength = 0.0CEP = 0.0}field-sampling{ sample = falsetype = at-pointfield = Exfield = Eyfield = Ezdirectory = ./base-name = field-sampling/fieldrhythm = 3.2position = (0.0, 0.0, 110.0)}field-visualization{
HAPTER 4. EXAMPLES 47 sample = falsefield = Exfield = Eyfield = Ezdirectory = ./base-name = field-visualization/fieldrhythm = 8.0}field-profile{ sample = falsefield = Qdirectory = ./base-name = field-profile/fieldrhythm = 100}}UNDULATOR{ undulator-type = opticalbeam-type = plane-waveposition = ( 0.0, 0.0, 0.0 )direction = ( 0.0, 0.0,-1.0 )polarization = ( 0.0, 1.0, 0.0 )strength-parameter = 0.5signal-type = flat-topwavelength = 1.0e3variance = 1200.0e3offset = 600082.0CEP = 0.0}FEL-OUTPUT{ radiation-power{ sample = truetype = at-pointdirectory = ./base-name = power-sampling/powerdistance-from-bunch = 82normalized-frequency = 1.000normalized-frequency = 2.000normalized-frequency = 3.000}}
Fig. 4.7a illustrates the radiation field 82 nm away from the bunch center. In addition, Fig. 4.7b shows the radiated power in termsof travelled undulator distance computed using MITHRA, illustrating the effect of space charge and energy spread. It is observedthat the gain obtained in this regime is very small compared with a usual static undulators. The reason for this effect is the very largeshot noise in the bunch because of the low number of particles in each micro-bunch. Note that in this simulation, each electron ismodeled as one single particle. The strong shot noise causes a strong initial radiation, which reaches the expected saturation powerafter a low gain. As a matter of fact, the micro-bunching process increases the coherence of the output radiation rather than poweramplification. Another aspect in this regime of interaction is the generation of strong higher order harmonics, which are depicted upto the third harmonic in Fig. 4.7c. Note that the accuracy of the results decreases for higher harmonics due to the required resolutionin the computational mesh. To show that the micro-bunching effect takes place in this regime as well, the bunching factor of the
HAPTER 4. EXAMPLES 48 z ( µ m) -3-2-10123 E y ( G V / m ) z ( µ m) l og P [ W ] without space-charge and energy spreadwith space-charge without energy spreadwith space-charge and energy spread ( a ) ( b ) z ( µ m) l og P [ W ] first harmonicsecond harmonicthird harmonic z ( µ m) B un c h i ng f a c t o r ( c ) ( d ) Figure 4.7: (a) Electric field of the generated radiation in front of the bunch, (b) the total radiated power measured at 82 nm distancefrom the bunch center in terms of the traveled distance, (c) the same radiation power for various harmonic orders, and (d) bunchingfactor of the considered bunch during the ICS interaction.electron beam is depicted in Fig. 4.7d. The bunching of the electrons due to the ICS interaction is clearly observed in the plotof bunching factor. According to the depicted power and pulse shape, total number of emitted photons is approximately equal to . × . To demonstrate the presented hypothesis related to the micro-bunching of bunches with low number of electrons perwavelength bucket, we perform an unreal simulation, where each electron is presented by 1000 particles. In this case, each particlerepresents a charge 1000 times smaller than the charge of one electron. In addition, we assume an initial bunching factor equal to0.001 for the input bunch to trigger the FEL gain. In Fig. 4.8, the radiation of such a charge configuration is depicted. The resultsclearly reveal the radiation start from much lower powers, possibility of achieving the FEL gain and saturating in the same powerlevel as above, thereby confirming the above theory for radiation of low density electron bunches. HAPTER 4. EXAMPLES 49 z ( µ m) -2-101234 l og P [ W ] Figure 4.8: The total radiated power measured at 82 nm distance from the bunch center in terms of the traveled distance for animaginary bunch where each electron is represented by a cloud of 1000 particles. ibliography [1] Peter Schm¨user, Martin Dohlus, and J¨org Rossbach.
Ultraviolet and soft X-ray free-electron lasers: introduction to physicalprinciples, experimental results, technological challenges , volume 229. Springer, 2008.[2] Evgeny L Saldin, E Evgeny A Schneidmiller, and Mikhail V Yurkov.
The physics of free electron lasers . Springer, 2000.[3] Henry P Freund.
Principles of free-electron lasers . Springer Science & Business Media, 2012.[4] A.E. Siegman.
Lasers . University Science Books, 1986. ISBN 9780935702118. URL https://books.google.de/books?id=1BZVwUZLTkAC .[5] AS Gilmour Jr. Microwave tubes.
Dedham, MA, Artech House, 1986, 502 p. , 1, 1986.[6] Eberhard Jaeschke, Shaukat Khan, Jochen R Schneider, and Jerome B Hastings.
Synchrotron Light Sources and Free-ElectronLasers . Springer, 2015.[7] S Reiche. Genesis 1.3: a fully 3d time-dependent fel simulation code.
Nuclear Instruments and Methods in Physics ResearchSection A: Accelerators, Spectrometers, Detectors and Associated Equipment , 429(1):243–248, 1999.[8] Sandra G Biedron, Henry P Freund, and Stephen V Milton. 3d fel code for the simulation of a high-gain harmonic generationexperiment. In
Optoelectronics’ 99-Integrated Optoelectronic Devices , pages 96–108. International Society for Optics andPhotonics, 1999.[9] Trach Minh Tran and JS Wurtele. TDA - a three-dimensional axisymmetric code for free-electron-laser (fel) simulation.
Computer Physics Communications , 54(2):263–272, 1989.[10] B Faatz, W Fawley, P Pierini, S Reiche, G Travish, D Whittum, and J Wurtele. Tda3d: Updates and improvements to thewidely used three-dimensional free electron laser simulation.
Nuclear Instruments and Methods in Physics Research SectionA: Accelerators, Spectrometers, Detectors and Associated Equipment , 393(1):277–279, 1997.[11] William M Fawley. A user manual for ginger and its post-processor xplotgin.
Lawrence Berkeley National Laboratory , 2002.[12] L Giannessi. Overview of perseo, a system for simulating fel dynamics in mathcad. In
Proceedings of the Free-Electron LaserConference , 2006.[13] A Bacci, C Maroli, V Petrillo, AR Rossi, L Serafini, and P Tomassini. Compact x-ray free-electron laser based on an opticalundulator.
Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors andAssociated Equipment , 587(2):388–397, 2008.[14] Roger J Dejus, Oleg A Shevchenko, and Nikolai A Vinokurov. An integral equation based computer code for high-gain free-electron lasers.
Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors andAssociated Equipment , 429(1):225–228, 1999.[15] EL Saldin, EA Schneidmiller, and MV Yurkov. Fast: a three-dimensional time-dependent fel simulation code.
NuclearInstruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment ,429(1):233–237, 1999.[16] Igor A Andriyash, Remi Lehe, and Victor Malka. A spectral unaveraged algorithm for free electron laser simulations.
Journalof Computational Physics , 282:397–409, 2015. 50
IBLIOGRAPHY 51 [17] LT Campbell and BWJ McNeil. Puffin: A three dimensional, unaveraged free electron laser simulation code.
Physics ofPlasmas (1994-present) , 19(9):093119, 2012.[18] SG Biedron, YC Chae, Roger J Dejus, B Faatz, HP Freund, SV Milton, H-D Nuhn, and S Reiche. Multi-dimensional free-electron laser simulation codes: a comparison study.
Nuclear Instruments and Methods in Physics Research Section A: Accel-erators, Spectrometers, Detectors and Associated Equipment , 445(1):110–115, 2000.[19] SPD Mangles, CD Murphy, Z Najmudin, AGR Thomas, JL Collier, AE Dangor, EJ Divall, PS Foster, JG Gallacher, CJ Hooker,et al. Monoenergetic beams of relativistic electrons from intense laser–plasma interactions.
Nature , 431(7008):535–538, 2004.[20] J´erˆome Faure, Yannick Glinec, A Pukhov, S Kiselev, S Gordienko, E Lefebvre, J-P Rousseau, F Burgy, and Victor Malka. Alaser–plasma accelerator producing monoenergetic electron beams.
Nature , 431(7008):541–544, 2004.[21] CGR Geddes, Cs Toth, J Van Tilborg, E Esarey, CB Schroeder, D Bruhwiler, C Nieter, J Cary, and WP Leemans. High-qualityelectron beams from a laser wakefield accelerator using plasma-channel guiding.
Nature , 431(7008):538–541, 2004.[22] T Tajima and JM Dawson. Laser electron accelerator.
Physical Review Letters , 43(4):267, 1979.[23] O Lundh, J Lim, C Rechatin, L Ammoura, A Ben-Ismail, X Davoine, Guilhem Gallot, Jean-Philippe Goddet, E Lefebvre,Victor Malka, et al. Few femtosecond, few kiloampere electron bunch produced by a laser-plasma accelerator.
Nature Physics ,7(3):219–222, 2011.[24] R Joel England, Robert J Noble, Karl Bane, David H Dowell, Cho-Kuen Ng, James E Spencer, Sami Tantawi, Ziran Wu,Robert L Byer, Edgar Peralta, et al. Dielectric laser accelerators.
Reviews of Modern Physics , 86(4):1337, 2014.[25] Emilio A Nanni, Wenqian R Huang, Kyung-Han Hong, Koustuban Ravi, Arya Fallahi, Gustavo Moriena, RJ Dwayne Miller,and Franz X K¨artner. Terahertz-driven linear electron acceleration.
Nature communications , 6, 2015.[26] Arya Fallahi, Moein Fakhari, Alireza Yahaghi, Miguel Arrieta, et al. Short electron bunch generation using single-cycleultrafast electron guns. arXiv preprint arXiv:1606.02153 , 2016.[27] FX K¨artner, F Ahr, A-L Calendron, H C¸ ankaya, S Carbajo, G Chang, G Cirmi, K D¨orner, U Dorda, A Fallahi, et al. Axsis:Exploring the frontiers in attosecond x-ray science, imaging and spectroscopy.
Nuclear Instruments and Methods in PhysicsResearch Section A: Accelerators, Spectrometers, Detectors and Associated Equipment , 2016.[28] Allen Taflove and Susan C Hagness.
Computational electrodynamics . Artech house publishers, 2000.[29] Peicheng Yu, Xinlu Xu, Viktor K Decyk, Weiming An, Jorge Vieira, Frank S Tsung, Ricardo A Fonseca, Wei Lu, Luis O Silva,and Warren B Mori. Modeling of laser wakefield acceleration in lorentz boosted frame using em-pic code with spectral solver.
Journal of Computational Physics , 266:124–138, 2014.[30] Jean-Luc Vay, Irving Haber, and Brendan B Godfrey. A domain decomposition method for pseudo-spectral electromagneticsimulations of plasmas.
Journal of Computational Physics , 243:260–268, 2013.[31] Jean-Luc Vay, DP Grote, RH Cohen, and Alex Friedman. Novel methods in the particle-in-cell accelerator code-frameworkwarp.
Computational Science & Discovery , 5(1):014019, 2012.[32] Kurt L Shlager and John B Schneider. Comparison of the dispersion properties of several low-dispersion finite-differencetime-domain algorithms.
IEEE Transactions on Antennas and Propagation , 51(3):642–653, 2003.[33] Bezalel Finkelstein and Raphael Kastner. Finite difference time domain dispersion reduction schemes.
Journal of Computa-tional Physics , 221(1):422–438, 2007.[34] Jay P Boris. Acceleration calculation from a scalar potential. Technical report, Princeton Univ., NJ Plasma Physics Lab., 1970.[35] JP Boris. Relativistic plasma simulation-optimization of a hybrid code. In
Proc. Fourth Conf. Num. Sim. Plasmas, Naval Res.Lab, Wash. DC , pages 3–67, 1970.[36] John David Jackson and Ronald F Fox. Classical electrodynamics.
American Journal of Physics , 67:841, 1999.[37] D Sagan, JA Crittenden, D Rubin, and E Forest. A magnetic field model for wigglers and undulators. In
Particle AcceleratorConference, 2003. PAC 2003. Proceedings of the , volume 2, pages 1023–1025. IEEE, 2003.
IBLIOGRAPHY 52 [38] C Pellegrini, A Marinelli, and S Reiche. The physics of x-ray free-electron lasers.
Reviews of Modern Physics , 88(1):015006,2016.[39] Sven Reiche. Numerical studies for a single pass high gain free-electron laser. Technical report, DESY, 2000.[40] C. Penman and B.W.J. McNeil. Simulation of input electron noise in the free-electron laser.
Optics Communications , 90(1):82– 84, 1992. ISSN 0030-4018.[41] C Maroli and V Petrillo. Effects of the low-frequency backward wave in high-gain free-electron lasers.
Optics communications ,183(1):139–147, 2000.[42] Trach-Minh Tran and JS Wurtele. Review of free-electron-laser (fel) simulation techniques. Technical report, 1990.[43] Eric Esarey, Sally K Ride, and Phillip Sprangle. Nonlinear thomson scattering of intense laser pulses from beams and plasmas.