The general particle tracer code : design, implementation and application

(1)

(2)

(3)

The General Particle Tracer code

Design, implementation and application

PROEFSCHRIFT

TER VERKRIJGING VAN DE GRAAD VAN DOCTOR AAN DE TECHNISCHE UNIVERSITEIT EINDHOVEN, OP GEZAG VAN DE RECTOR MAGNIFICUS, PROF.DR. M. REM, VOOR EEN COMMISSIE AANGEWEZEN DOOR HET COLLEGE VOOR PROMOTIES IN HET OPENBAAR TE VERDEDIGEN OP

MAANDAG 2 APRIL 2001 OM 15.00 UUR

DOOR

BAS VAN DER GEER, GEBOREN TE SOEST EN

(4)

prof.dr. M.J. van der Wiel en

prof.dr.ir. T.J. Schep Copromotor: dr. J.I.M. Botman

CIP-DATA LIBRARY TECHNISCHE UNIVERSITEIT EINDHOVEN van der Geer, S.B.

de Loos, M.J.

The General Particle Tracer code : Design, Implementation and Application / by S.B. van der Geer, M.J. de Loos. –

Eindhoven: Technische Universiteit Eindhoven, 2001. – Proefschrift. – ISBN 90-386-1739-9

NUGI 812

Subject headings: particle optics / software engineering / electron beams / free electron lasers

(5)

(6)

(7)

v

2.10.6 Broydens method 59 2.10.7 The optimizer 59 3 GPT code design 63 3.1 Introduction 63 3.2 General considerations 65 3.2.1 Efficiency 65 3.2.2 Flexibility 65 3.2.3 Scalability 66 3.2.4 Diagnostic messages 67 3.2.5 Multi-platform 67 3.2.6 Reliability 68 3.2.7 Convenience 68 3.2.8 Working in a team 69 3.3 Language 69 3.4 Inputfile 70 3.5 Custom elements 71

3.5.1 The initialization routine 72

3.5.2 Custom element interface 74

3.5.3 Callback functions 75

3.5.4 The info structure 76

3.5.5 Reading parameters into the info structure 78

3.5.6 Calculating electromagnetic fields 78

3.5.7 ODE advanced callback functions 80

3.5.8 Removing a particle 83

3.5.9 Multiprocessing on the FPR interface 83

(9)

vii

3.6 GDF 85

3.6.1 Introduction 85

3.6.2 Disk format and driver programs 87

3.6.3 The GDF library and its memory format 89

3.6.4 Conversion programs 91

3.6.5 GDFA data analysis 92

3.7 GPTwin 93

3.7.1 MFC 95

3.7.2 Running GPT 95

3.7.3 Plotting 96

3.7.4 Elements and compilation 98

4 Design of a 100 fs photo-gun 101 4.1 Introduction 101 4.2 Design process 102 4.3 Diode 104 4.3.1 Required field 104 4.3.2 Diode set-up 106 4.3.3 GPT simulations 108

4.3.4 Reference simulation results 120

4.3.5 Electrostatic compensation of non-linear space-charge effects 124

4.3.6 Initial beam radius 128

4.3.7 Focusing 131

4.3.8 Laser parameters 133

4.3.9 Scan of beam charge 136

4.3.10 Conclusion 137

4.4 The rf booster 138

4.4.1 Required accelerator field 139

4.4.2 The modified BNL design 140

4.4.3 Superfish calculations 142

4.4.4 Resonant frequencies 145

4.4.5 RF incoupling 148

4.5 Combined diode and rf booster simulation results 153

4.5.1 Set-up 153

4.5.2 Optimization 156

4.5.3 Short bunch, low emittance, divergent beam 157

4.5.4 Parallel beam 160

4.5.5 Waist after accelerator, emittance compensation 161

(10)

5 The energy recovery system of the ‘Rijnhuizen’ Free Electron Maser 165

5.1 Introduction 165

5.1.1 Beam line 165

5.1.2 Efficiency 167

5.1.3 Simulations 169

5.2 FEL interaction in the undulators 170

5.3 Beam transport downstream of the undulator 173

5.3.1 High-energy transport section 174

5.3.2 Decelerator 175

5.3.3 Low-energy transport section 176

5.3.4 Simulations 181

5.4 The depressed collector 183

5.4.1 The collector 184

5.4.2 Scattering 187

5.4.3 The magnetic deflection system 190

5.4.4 Current dissipation and return current 195

5.4.5 Power dissipation 197

5.5 Conclusion 201

References 203

Publications related to this thesis 207

Summary 209

Samenvatting 211

Curriculum vitarum 213

Dankwoord 215

(11)

1 Introduction

The high-energy physics community, users of the most advanced accelerators, numbers about ten thousand scientists and engineers. Over 10 billion USD is spent every year on medical linear accelerators alone. Synchrotron radiation sources produce more than a million user-hours of beam time per year. Clearly, charged particle beams are ‘big business’ and essential tools for scientific, medical and industrial use.

Numerical simulations are essential for the design and understanding of charged particle accelerators. This thesis describes the design and application of a new code for charged particle beam simulations, developed by us within our own company Pulsar Physics. Apart from implementing and marketing this code commercially, we have used the code for the design of two novel and challenging electron beam experiments described in this thesis:

• A photo-gun delivering 100 fs high-brightness electron bunches • A beam and energy recovery system in a Free Electron Maser

1.1 Charged particle beam simulations

Devices producing charged particle beams consist of a source emitting the particles, followed by an accelerator. The principles are relatively simple: To emit electrons for example, it is sufficient to shine light on or heat a metal surface. Protons and ions can be extracted from ionized gas. The accelerator itself can consist of just two plates with different electric potential. Although all these techniques are still commonly in use today, science, industry and the medical community demand a very high level of sophistication. As many particles as possible should be emitted, accelerated to higher and higher energies, focused to tiny spot sizes with extra long or unimaginably short pulses. And, of course, this should all be not too expensive.

Real-life machines often consist of one or more accelerating sections, sometimes circular, optionally ending in a sophisticated transport section to deliver the beam with high accuracy to the target location. The design of these machines relies heavily on predicting the beam behavior by analytical expressions and numerical simulations. The basic equations of motion for charged particles in electromagnetic fields are well known and relatively easy to solve. The interaction of the particle beam with itself due its own charge, known as space-charge, is also based on relatively simple equations. Unfortunately however, the space-charge equations are hard to solve because every

(12)

particle, and there can be billions of those, interacts with every other particle. And to make things even worse, the interaction of the particle beam with surrounding material such as beam pipes or accelerating sections requires the calculation of the interaction between every particle and every individual component in the set-up.

Due to the space-charge and the interaction with the surrounding material, it is virtually impossible to perform a complete numerical simulation including all physical effects present in an accelerator. Therefore, simplifications need to be made. A number of common approaches made in popular general-purpose codes can be distinguished: • The most radical approximation assumes all transverse and longitudinal forces to

be linear with the distance from the trajectory of a reference particle. Every beam line component is written as a matrix mapping the beam parameters before and after the element. The TRACE3D [1] code is a relatively modern example capable of handling space-charge forces linear with transverse distance. Apart from initial designs, this approximation is in many cases a too simplistic view.

• The linear approach can be generalized to higher order expansions. The TRANSPORT [2] and MAD [3] codes for example can be run up to 3rd order and the COSY [4] code can be used to any desired order. The slight loss of computation speed is compensated by far more accurate results compared to the linear case. For circular machines and many beam line designs these codes are very useful. However, they are not capable of modeling the space-charge effects in an arbitrarily shaped charged particle bunch.

• A number of sample particles representing the complete bunch are tracked the in time domain through the electromagnetic fields of the set-up without any expansion into orders. Space-charge effects taking into account the shape of the particle beam can either be calculated using particle-particle interaction or using a mesh mechanism. Unfortunately, the particle-particle method introduces statistical noise and the mesh method must make some unrealistic assumptions. Still, the tracking algorithms are much more powerful than the matrix formalism at the cost of a significant increase in CPU time. The PARMELA [5] code is currently the most widely used program in this category.

• Strong interaction with walls is a typical situation at the particle source, the gun. Solving the electrostatic potential of the beam and the gun iteratively while tracing a sample beam until a stable solution is found is the approach used in the two dimensional EGUN [6] code. This approach only works for continuous beams or relatively long bunches.

• Time dependent interaction with walls is very time-consuming to solve numerically. The three-dimensional MAFIA [7] code for example subdivides both the beam and the accelerator in a 3D mesh and solves all interactions self-consistently while tracing sample particle trajectories. Almost all relevant physical effects can be included but, because the mesh must be quite dense to avoid

(13)

GPT: A new beam line design code 3

unphysical effects, the user’s patience will be tested and the method can only be applied to relatively short structures.

Because of the lack of a fast, all-purpose beam line design code, many codes are written using ‘intermediate’ approaches for specific situations combining numerical methods with analytical solutions. The HOMDYN code for example slices a bunch into a number of longitudinal slabs to be able to calculate so-called emittance compensation schemes while the speed is maintained using analytical expressions based on the envelope of fixed distributions for the cylindrically symmetric transverse part. Such custom-codes are essential to test new algorithms, but not suited for generic accelerator design due to the approximations made.

1.2 GPT: A new beam line design code

Because of the enormous advances in computer technology, it is often possible to perform more accurate beam line simulations than attainable using the matrix formalism. On the other hand, true time dependent and 3D electromagnetic design includes all wall interactions and is often not necessary and too time consuming. In many design situations, the particle tracking approach is therefore a good compromise between the accuracy of the results and the CPU time spent. For these reasons, a large number of scientists use particle tracking techniques as their main design tool for accelerators and beam lines for a wide variety of applications.

The requirements for a good particle tracking code are quite general. Unfortunately, the only general purpose code available in the 1990’s was PARMELA. Although this package is in many aspects still up-to-the-job it has been called a dinosaur for good reasons. When the code was developed a few decades ago it made the best use of the computing and language power available of that time, but in the early 1990’s, the time seemed right for a remake.

PARMELA simulations for the design of the injector of the free electron laser FELIX [8] at the FOM Institute for Plasma Physics ‘Rijnhuizen’ where both frustrating and time-consuming. Because our attempts to modernize PARMELA were unsuccessful, we started the General Particle Tracer (GPT) project in 1992. We anticipated a worldwide demand for something better than PARMELA, founded the company Pulsar Physics [9] and started developing the GPT code as a commercial product.

Similar to PARMELA, the GPT code is based on 3D tracking of charged particle trajectories, but GPT uses a higher-order solver with automatic stepsize control. Arbitrary positioning in 3D of all beamline components, the ability to add new

(14)

space-charge models, better output options and the use of SI units throughout the code were just a few of the additional wishes for this new code. To make the GPT code useful to the scientific public, it was anticipated that it would be impossible for us to create every element users might require in the future. Therefore, an object-oriented approach was created to allow custom beam-line components, specifying the electromagnetic fields as function of position and time, to be developed easily and to be exchanged with other users.

Table 1-A: Institutes worldwide with a commercial GPT license. Not listed in the table are GPT versions used in different departments within the same institute and updates.

Institute Name Country Year

Stanford University USA 1994

FOM-Rijnhuizen Netherlands 1995

University of Abertay Dundee UK 1996

Sumitomo Heavy Industries, Ltd. Japan 1996

University of Strathclyde UK 1996

Fermi National Accelerator Laboratory (FNAL) USA 1996

Forschungszentrum Rossendorf (FZR) Germany 1997

Pohang Accelerator Laboratory (POSTECH) South Korea 1997

Technische Universität Darmstadt Germany 1997

Korea Atomic Energy Research Institute (KAERI) South Korea 1997

Japan Atomic Energy Research Institute (JAERI) Japan 1997

Deutsches Electronen-Synchrotron (DESY) Germany 1997

Groningen University Netherlands 1997

Eindhoven University of Technology (TUE) Netherlands 1998

Rafael Laboratories Israel 1998

Ishikawajima-Harima Heavy Industries Co., Ltd. Japan 1998

Los Alamos National Laboratory (LANL) USA 1998

Tel Aviv University Israel 1999

Océ Printing Systems Germany 2000

European Synchrotron Radiation Facility (ESRF) France 2000

Ankara University Turkey 2000

Utrecht University Netherlands 2000

Marshall Space Flight Center (NASA) USA 2001

Only two years after the start of the GPT project, the first commercial user from Stanford University received the first GPT version in 1994. All initial design goals, both the computational physics and the software engineering aspects as described in chapters two and three of this thesis respectively, had been achieved. Unfortunately for the FELIX injector design, the first version of the GPT code came too late to be of much use. However, as can be seen in Table 1-A, many institutes around the world

(15)

Towards a table-top XUV source 5

have used GPT to solve their problems. Later additions to the GPT code include free-electron laser interaction, more space-charge models, 3D scattering, more built-in elements and various conversion utilities to- and from other codes.

Until today we have developed about two new versions of GPT each year, but our attention has somewhat shifted to using, and adapting where necessary, GPT for various projects. These projects, including the ones described in chapters four and five in this thesis have all been carried out on contract basis with our company Pulsar Physics.

1.3 Towards a table-top XUV source

At the Eindhoven University of Technology (TUE), a project has been started aiming at the production of ultra-short radiation pulses generated by relativistic electrons. The long term goal of the project is to create a table top free-electron laser for (X)UV-radiation, producing intense 100 fs pulses that can subsequently be used in various disciplines of scientific research. The first step is to achieve 100 fs short electron bunches. Our contribution to the project, as described in chapter four of this thesis, is the electromagnetic design and electron beam trajectory calculations with GPT in the first two stages of the set-up: DC acceleration in a 1 GV/m field followed by a state-of-the-art high-gradient radio-frequency (rf) acceleration section. At the time of writing this thesis, the first experimental tests of the set-up are being performed.

Synchrotron radiation facilities throughout the world combined produce over a million hours beam time. The past three generations of these ‘science-factories’ have increased peak intensity by many orders of magnitude, but still the demand for even higher intensities and shorter pulses is stronger than ever. High intensity, 100 fs light pulses, combined with short wavelengths around 0.1 nm, or 1 Å, opens the field of stroboscopic images of vibrating atoms. It would allow researchers to ‘see’ the dynamics of molecules. However, as the synchrotron user facilities now operate very close to their theoretical limits, their evolution seems to be coming to an end and a different approach is required.

The most promising scheme is a linear electron accelerator (linac) followed by a Self-Amplified Spontaneous Emission (SASE) Free Electron Laser (FEL). The first device working on this principle is the DESY project in Germany. In this set-up, a relatively long electron bunch of 2 ps is generated by photo-emission from a semiconductor surface. This bunch is accelerated by roughly two hundred meters of superconducting rf structures to sufficiently high energy. Properly compressing the bunch to the required 100 fs will be a complicated process due to emitted Coherent Synchrotron Radiation

(16)

(CSR), degrading the bunch quality. After acceleration and compression, the beam is sent into a periodic alternating magnetic field, an undulator. In this undulator, with a length of tens of meters, part of the energy of the electron bunch is converted to VUV-radiation of about 6 nanometer with a single but selectable wavelength. Because in the VUV and X-Ray wavelength range there are no mirrors that can be used to send the light back and forth as in a conventional FEL, the radiation must be generated in a single pass within the undulator. The lack of a suitable ‘drive-laser’ makes starting from noise, SASE, the only option, hereby further increasing the requirements of the machine.

The DESY machine is very large, very expensive and an X-ray version operating at 0.1 nanometer even will be worse. To develop a compact alternative, i.e. to miniaturize the accelerator, much higher acceleration gradients are needed. Evolution in current acceleration techniques, better superconducting technology and higher rf frequencies could reduce the accelerator size, but to reduce it to a few meters requires a radically different approach. Theoretical predictions and first experimental results point towards plasma acceleration schemes. A promising option is laser wakefield acceleration, where acceleration gradients of over 100 GV/m are achievable theoretically.

Because most plasma acceleration schemes require the input of pre-accelerated, ultra-short, high-brightness electron bunches, developing an electron source for bunches of 100 fs is a good starting point for future research on compact radiation devices. In the proposed TUE set-up a 100 fs electron bunch is first accelerated in a quasi-DC, 1 ns, 1 GV/m acceleration gap. This accelerates the bunch before it explodes due to its own charge and eliminates the need for downstream compression, avoiding the CSR problems. Although the long-term goal of the TUE project is a table top (X)UV laser, the 100 fs photo-gun followed by a booster rf accelerator can already be used for interesting experiments such as the production of 100 fs transition or Cherenkov radiation or Compton scattering on a laser pulse to produce hard X-rays. The design and optimization process for this TUE short-bunch, low-emittance photo-injector is presented in chapter four.

1.4 Energy recovery in a free electron maser

Free Electron Lasers (FELs) are a reliable source of coherent, high-power and continuously tunable radiation. As mentioned before, they produce this radiation by forcing relativistic electrons into a wiggling motion in an undulator. Unfortunately, because only a few percent of the energy of the relativistic electrons is converted into radiation using this mechanism, the efficiency of FELs is typically not higher than that of ‘regular’ lasers. However, by decelerating and collecting the spent electron beam,

(17)

Energy recovery in a free electron maser 7

and thus recovering its energy, high efficiency on the order of 50% can be added to the impressive list of characteristics of FEL radiation.

The FOM-Fusion Free Electron Maser (FEM) is an FEL-like device currently under construction at the FOM-Institute for Plasma Physics ‘Rijnhuizen’ in the Netherlands. Although the FEM device is scientifically very interesting, the main drive for the development is its application for future fusion devices for heating and current drive at the electron-cyclotron frequency in magnetically confined plasmas such as ITER. Compared to gyrotrons, which have comparable efficiencies, the FEM has the advantage of fast broadband tunability, which is essential for e.g. the control of instabilities in a tokamak plasma. Other applications involve power beaming to communication satellites and certain military uses.

For the first time ever, the FEM will demonstrate high-efficiency at power levels around 1 MW during relatively long pulses of 100 ms. Our contribution to the FEM is the design of the beam- and energy recovery system, one of the critical parts of the machine. This system consists of a transport section and a multi-stage collector where the electrons are absorbed on the backside of several collector plates. Although the concept of such a collector is not new, the design requirements for this collector with its high power are extreme in the sense that even one percent return current in the FEM or a non-uniform electron distribution on the collection plates will severely damage the machine. Where typical collection systems are cylindrically symmetric, we designed a 3D off-axis bending scheme to reduce the return current to the required level while improving the uniformity of the power and current distribution on the collector plates. Our design, presented in chapter five, heavily relies on GPT simulations predicting the electron beam behavior in the complete set-up. Interaction between the electrons and the radiation wave, space-charge, high-current low-energy beam transport and 3D scattering inside the collector are all included. To model these processes, two new extensions have been designed and implemented for the GPT code: Particle-wave interaction and scattering of particles off surfaces. Currently both extensions are commercially licensed to different institutes for similar design work. Until today we are not aware of any other code capable of simulating electron beam behavior in the complete FEM, including the 3D collector.

(18)

(19)

2 GPT: Physics and

Mathematics

2.1 Introduction

The General Particle Tracer (GPT) is a software package developed to study 3D charged particle dynamics in electromagnetic fields. Developed as a commercial product, GPT is being used worldwide by a large number of scientific institutes, mainly for accelerator and beam line design. This chapter describes the ‘equations and numerical algorithms of GPT’. The software-engineering details of GPT are presented in chapter three.

GPT tracks any number of charged macro/sample-particles through electromagnetic fields, taking all 3D effects and space-charge forces into account. Due to the general character of the code and its flexibility, GPT is suited for a large number of applications. To explain the GPT code, it will mainly be compared with the world standard PARMELA [5] code, developed at Los Alamos National Laboratory. Despite significant differences in underlying physics and implementation, from a user point of view the GPT and PARMELA codes are comparable.

GPT solves the equations of motion for the macro-particles relativistically in the time-domain using a 5th order embedded Runge-Kutta integrator with adaptive stepsize control. This allows the user to select the required accuracy, while the simulation time is always kept to a minimum. PARMELA uses a lower order scheme without monitoring its internal accuracy. GPT allows the equations of motion to be combined with additional differential equations to solve for example beam-loading in a linac or FEL interaction self-consistently.

Both GPT and PARMELA come with a large set of standard elements for basic structures such as solenoids, quadrupoles and accelerating structures. Interfaces with external Poisson and rf solvers calculating electromagnetic fields in various geometries are provided using field-maps. Custom elements or fringe-fields for specific cases can easily be added and elements can be positioned anywhere and in any direction in 3D space.

(20)

The self-fields of the particles are also part of the electromagnetic fields through which the particles are tracked. In GPT, space-charge effects can be calculated using various models depending on the type of simulation. Cylindrically symmetric or (semi)continuous beams are best handled with the 2D space-charge models. They are based on relativistic point-to-ray and point-to-circle interaction. The 3D routine consists of a relativistic point-to-point model. PARMELA has a 2D mesh based routine for cylindrically symmetric beams and a 3D point-to-point model.

GPT has two available output modes: time and position output. Time output writes all particle coordinates at user specified time(s). Position output writes all particle coordinates passing any plane in 3D space. This output mode is also known as ‘non-destructive screen’ output. Optionally, the electric and magnetic fields at the particle coordinates are also output. PARMELA has end-of-element output only and does not provide the electric and magnetic field for diagnostic purposes, making it impossible to investigate particle dynamics inside a beam line component.

Being able to simulate charged particle dynamics in time-dependent 3D electromagnetic field configurations is usually far from sufficient for serious accelerator and beam line design. The simulation data must be analyzed, parameters must be scanned and typically a comparison must be made between different scenarios. GPT is accompanied by a number of pre- and postprocessing programs, a scanning utility, a multi-dimensional solver and constraint optimizer and various interfaces to other software packages to ease this design process. These tools are all combined with GPT and integrated into the Microsoft Windows based graphical user interface, GPTwin.

2.2 Equations of motion

2.2.1 Introduction

GPT calculates the trajectories of the charged particles through the combined electromagnetic fields of the set-up and the self-fields of the beam in time domain. This is the same principle as used in PARMELA. Because it is not possible to solve these trajectories for all elementary particles in a typical beam, both GPT and PARMELA group particles together to form macro-particles. Because the mass-to-charge ratio of these macro-particles is identical to the elementary particles, the equations of motion are identical. Depending on the application, a few hundred to a few thousand macro-particles are sufficient to represent the beam while over a billion elementary macro-particles is not uncommon.

(21)

Equations of motion 11

The relativistic equations of motion for (macro) particle i are given by:

(

)

1 2₊ = = = × + = i i i i i i i i i c c dt d mc q dt d β β β v x B v E β γ γ γ [2-1]

where the position x and the normalized momentum γβ=p/mc are used as the coordinates of a particle. All equations of motion are solved in the laboratory frame. Various methods to calculate the E and B fields at the position of the particles are described in sections 2.4 and 2.5.

A difference between PARMELA and GPT is the fact that GPT uses SI units. As a result, time is specified in seconds using GPT, not in degrees of the phase of an oscillator as in PARMELA. Electric and magnetic fields are specified in volt per meter and Tesla respectively, as opposed to volt per centimeter and gauss.

The equations of motion can not be solved for each particle individually because, due to the space-charge, the force on each particle depends on the position of all other particles. Therefore we introduce a notational vector y(t) containing the six coordinates of all the particles as function of time t. The equations of motion [2-1] can then be rewritten as a simple Ordinary Differential Equation (ODE):

(

)

d t dt t t y f y ( ) , ( ) = [2-2]

where f

₍

t, ( )yt

₎

contains the equations [2-1] for every particle. The boundary

conditions are the particle coordinates at a specified time.

A number of standard textbook methods exist for solving [2-2]. They are all based on taking finite steps in t while repeatedly evaluating f at various intermediate positions. Typically, the smaller the steps, the better the results at the cost of an increase in CPU time. Directly solving [2-2] however can be quite inefficient when (near) discontinuous fields are present at the beginning and end of beam line components. The following sections describe the PARMELA and GPT algorithms to solving [2-2].

2.2.2 The PARMELA approach

The main integration method of PARMELA is the leapfrog scheme, also known as drift-impulse-drift. It is the most simple integration method that guarantees a reasonable accuracy. The position and momentum variables are updated one after the other in time domain with a half step in between, resulting in second-order accuracy. In equations:

(22)

t t t f t c x x t dt d t i i i i i i ∆ ∆ + + = ∆ ⋅ + = = + + + + + ) , ( , | ) , ( 1 2 1 2 3 2 1 1 x β β β x β x f γ γ γ [2-3]

A historical advantage of the leapfrog scheme is the fact that it does not require additional storage. The main reason for using the method today is not high accuracy or efficiency, but stability. Because of the time symmetry, reversing time leads back to the starting point, the scheme is very insensitive to systematic errors.

Integrating from a to b is identical to integrating from a to some midpoint and from the midpoint to b. When the midpoint happens to be at the only discontinuity in the interval, the resulting two integrals are continuous and can be evaluated much more efficiently. PARMELA uses this approach by integrating to the boundaries of beam line components. Because element boundaries are defined in position space, not in time-domain, the required midpoint calculation can only be applied efficiently when the particle velocity is constant. In the leapfrog scheme this is possible, but the mechanism severely limits the choice on Ordinary Differential Equation (ODE) solvers to be used. Furthermore, true field-discontinuities do not exist in vacuum. When a beam line component is modeled properly, continuous fringe fields are always present. Because the discontinuity is typically introduced by using too simplistic models for beam line components, optimizing the ODE solver for these discontinuities does not solve the actual problem. Finally, space-charge effects bind the differential equations of all particles to each other. When the stepsize for one particle is reduced to the end of a beam line component, all other particles should make an identical step. This either results in the smallest stepsize dictating the overall step, or assuming the change in space-charge negligible in sub-step timescales.

Electromagnetic fields cause a change in momentum. However, when a particle is inside an element containing a magnetic field only, the total energy will not be affected. As a result, the equations of motion can be solved more accurately than solving [2-2] without this additional information. For this reason, PARMELA directly modifies γβz

to enforce conservation of energy in magnetostatic elements:

(

)

(

)

(

2 2 2

)

y y old, x x old, new y x old , , x y z x x z y z z y β β β β β β B B mc q β B B mc q β γ γ γ γ γ γ γ γ β β γ β β γ ∆ − ∆ − ∆ + ∆ + = − = ∆ − = ∆ β β [2-4]

(23)

As an example, PARMELA enforces the total momentum to be constant in a wiggler, an alternating magnetic structure. This increases the accuracy of the simulation, but the use is very limited because particles typically lose energy to the radiation field present inside a wiggler. As a result, PARMELA is not able to simulate a radiation field combined with a wiggler.

2.2.3 The GPT Runge-Kutta based ODE solver

Although discontinuities in electromagnetic fields are typically an over-simplification of the set-up, they have to be taken seriously because it is very convenient to work with hard-edge approximations in the early design stages. This eases interpretation of the results and allows the results to be compared with analytical/optical calculations. On the other hand, a higher order integration scheme is better suited for smooth fields. The fixed and user-defined stepsizes of PARMELA would ideally be changed into an algorithm that automatically chooses the correct stepsize. And in theory, the (near) discontinuities are also automatically taken care of when the stepsize is chosen automatically. We decided to follow this path and generically solve [2-2] without enforced boundary midpoints. It is the most general approach and it does not need to make any approximations about the space-charge fields. An additional advantage of generically solving [2-2] is the fact that additional variables and equations can easily be added to y(t) and f respectively. These extra differential equations can be a function of all the particle coordinates, time and the additional variables. Because [2-2] is solved simultaneously for all particles coordinates and extra variables, the results are always self-consistent. For example, beamloading in a linac and particle-wave interaction as described in sections 2.4.1 and 2.4.6 respectively, make use of this mechanism.

Because nobody likes to wait for a simulation code any longer than strictly necessary, the requirements for our generic ODE solver are clear:

• As few evaluations of f as possible.

• Discontinuities must be solved by adaptive stepsize control.

• Increasing the specified accuracy must converge to the analytically correct solution.

(24)

To test different integration routines we decided to run two simple test cases with the following textbook ODE integration routines:

Midpoint A half-step is taken using the derivative information at the start. The derivative information halfway is used to take a full step starting at the start. We implemented adaptive stepsize control by taking every step twice: Once as a full and once as two halve steps. The difference between these two is used as accuracy indication and improves the method from 2nd to 3rd-order.

Bulirsch-Stoer A number of midpoint sequences with gradually decreasing stepsizes is used to extrapolate the solution to zero stepsize. This extrapolation simultaneously produces an accuracy indication.

Runge-Kutta Six intermediate steps are taken to achieve a 5th_{-order estimate. The}

difference with an embedded 4th-order estimate is used as accuracy indication.

The first test is an integration of 10 periods of a sinus, similar to a particle trajectory inside a wiggler. The results are shown in Figure 2-1. To make the plots better readable, the imposed accuracy per step is given by 10–‘specified accuracy’_{. The achieved accuracy is}

defined as –10log(|result|), the number of correct digits in the final result. From the results it is clear that dynamic stepsizing works well for all algorithms because the achieved accuracy is always close to the specified accuracy. The main difference is the number of function evaluations as function of achieved accuracy. In this first test the Bulirsch-Stoer is superior by far especially in the high-accuracy range, as could be expected by the smoothness of the test function. Runge-Kutta and Midpoint require roughly two or three times more function evaluations for the same accuracy.

The second test, 10 discontinuous square blocks, produces very different results as shown in Figure 2-2. Bulirsch-Stoer is known to be an unwise choice for discontinuous functions. Without warning and with a reasonable number of function evaluations, the algorithm just comes up with the wrong result. Midpoint does not do any better. Runge-Kutta requires an order of magnitude more function evaluations compared to the sinus case, but it is the only one to succeed.

Clearly, the ‘best’ ODE solver cannot be defined. This strongly depends on the function to be integrated. For smooth functions and high accuracy Bulirsch-Stoer seems to be an ideal choice. However, because we need to be able to handle discontinuities and a precision of two or three digits is typically sufficient, we provided GPT with a fifth-order embedded Runge-Kutta ODE solver.

(25)

0 1 2 3 4 5 6 7 8

Midpoint specified accuracy

0 2 4 6 8 A chi ev ed a cc u ra cy 0 1 2 3 4 5 6 7 8

Midpoint achieved accuracy

0 500 1000 F unc tio n ev al u a tions 0 1 2 3 4 5 6 7 8

Bulirsch-Stoer specified accuracy

Bulirsch-Stoer achieved accuracy

Runge-Kutta specified accuracy

Runge-Kutta achieved accuracy

0 500 1000 F unc tio n ev al u a tions

Figure 2-1: Number of significant digits in the integration of 10 periods of a sinus using the Midpoint, Bulirsch-Stoer and Runge-Kutta integration algorithms.

(26)

0 1 2 3 4 5 6 7 8

Midpoint specified accuracy

Midpoint achieved accuracy

Bulirsch-Stoer specified accuracy

Bulirsch-Stoer achieved accuracy

Runge-Kutta specified accuracy

Runge-Kutta achieved accuracy

0 5000 10000 F unc tio n ev al u a tions

Figure 2-2: Number of significant digits in the integration of 10 blocks using the Midpoint, Bulirsch-Stoer and Runge-Kutta algorithms.

(27)

Coordinate systems 17

2.2.4 Additional differential equations

Apart from solving the individual particle trajectories, GPT has the capability to solve additional differential equations. This is not possible using PARMELA. The required additional variables are added directly to the vector y in [2-2]. Corresponding equations added to f in [2-2] can be a function of the additional variables as well as all particle coordinates. The results are always self-consistent with the particle trajectories because all equations are solved simultaneously.

The possibilities provided by this feature are endless, but the following two are further discussed in this thesis:

• Beamloading in a travelling wave linac as described in section 2.4.1. • Particle wave interaction in a HE waveguide mode as described in 2.4.6.

2.3 Coordinate

systems

PARMELA uses incremental positioning along the z-axis for the specification of the locations of the elements in the set-up. This assumes all elements to be side-by-side where drift spaces can be simulated using a special drift element without electromagnetic fields. Misaligned and off-axis positioned beam line components can be simulated, but 3D geometries as simple as two 90 degree rotated solenoids used as bending magnet can not be specified.

GPT uses absolute positioning to allow full control over misaligned and off-axis positioned elements. The user has the full six degrees of freedom in position and orientation. Drift spaces are not required, but can be inserted when clipping at a specified radius is desired.

The superposition principle allows the electromagnetic fields of all the elements in the set-up to be added when no interaction between elements is present. An all other cases, an external field solver must be used to calculate the correct combined field that can subsequently be imported as field-map.

As an optimization when many elements are present, GPT differentiates between two kinds of elements: Global elements and local elements. Global elements have a relatively long working range, i.e. a large solenoid, and all fields of all the global elements are added to obtain the total field. Local elements have a relatively short working range and their ranges are expected not to overlap. As a result, when a particle is in range of one local element, all the other local elements can be ignored to reduce CPU time.

(28)

2.3.1 Element Coordinate System (ECS)

GPT allows any standard beam line component to be arbitrarily positioned and oriented in 3D. To ease the development of the elements and to provide a consistent inputfile specification, the GPT kernel handles all necessary coordinate transforms while the element code can be written in any convenient coordinate system. For example, the code for a single-turn solenoid calculates its fields as function of current and radius in the most convenient coordinate system, centered on the origin with the z-axis as its normal. This private coordinate system, the Element Coordinate System (ECS) can be freely chosen by the developer for every element, but must always be a right-handed orthonormal Cartesian coordinate system. It is typically chosen to be in the center of the element aligned with the z-axis in the direction of the beam passage.

The base coordinate system of GPT, the World Coordinate System (WCS), is also an orthonormal right-handed Cartesian coordinate system. The relation between WCS and ECS is given by an orthonormal matrix M and an offset o as follows:

rWCS =MrECS+o [2-5]

or

(

)

rECS =M r−1 WCS−o [2-6]

where r is a coordinate measured in either WCS or ECS.

All particle positions are stored relative to WCS. To obtain the electromagnetic fields of every element, the particle coordinates are first transformed to the ECS of the element using [2-6]. Then, the fields are calculated in the ECS of the element, before they are transformed to WCS using:

E ME B MB WCS ECS WCS ECS = = [2-7]

2.3.2 Custom Coordinate System (CCS)

Determining the position and orientation of beam line components downstream a bending magnet is difficult in the WCS system. However, the location can easily be specified relative to a coordinate system with the z-axis aligned with the beam after the bend. For such situations, any number of additional orthonormal Custom Coordinate Systems (CCS) can be defined and used for the positioning of elements. For convenience, GPT output particle coordinates and fields can also be written relative to any CCS. The coordinate transformations relating the particle coordinates between WCS and a Custom Coordinate System, CCS, are similar to [2-5] and [2-6].

(29)

Physical models of selected elements 19

The ECS to WCS transforms are calculated by the GPT kernel by multiplying the transform from WCS to CCS and the transform from CCS to ECS:

M M M o M o o = = + CCS ECS CCS ECS CCS [2-8] They are used as in [2-5] and [2-6].

Because matrix-vector multiplications are costly in terms of CPU time, they are not performed for identity transformations. The code checks for these identity transformations during the initialization, thus avoiding the coordinate transformation overhead for straight beamline sections completely.

2.4 Physical models of selected elements

It is relatively simple to develop a custom GPT element when the electromagnetic field equations are known. For this reason, GPT has a large number of built-in elements and perhaps an even larger number of custom elements developed by the GPT user community. As a result, it is impossible to explain the physics of all included elements. This section describes a number of elements that are either particularly interesting or used in the projects described in chapters four and five.

2.4.1 Traveling wave accelerator with beamloading

The trwlinbm element models a constant gradient traveling wave buncher or linac with beamloading. The beamloading fields, caused by interaction between the particle beam and the accelerating structure, reduce the acceleration gradient as a direct consequence of Lenz’s Law. Although PARMELA is capable of calculating beamloading effects, the used assumptions only hold when the actual rf phase of the buncher or linac is close to the design phase. To investigate if the second linac of the FELIX accelerator [10,11] could be used as a decelerator, this assumption does not hold and the GPT element trwlinbm was developed. This element makes no assumptions about the phase of the particles with respect to the design phase by solving the phase and amplitude of the beamloading wave while tracing the macro-particles. The parameters of trwlinbm are listed in Table 2-A.

(30)

Table 2-A: Parameters of the trwlinbm element.

Parameter Description

ECS Element Coordinate System.

α0 Initial attenuation constant.

Rs Shunt impedance [Ω/m].

P0 Input power, design value [W].

P Input power, actual power used [W].

Ib Beam current [A].

γ0 Design gamma at entrance.

θ0 Bunch phase with respect to wave, design value.

ϕ RF phase offset.

ω Angular frequency [s-1_].

L Length of the structure [m].

An rf accelerating wave with constant amplitude travels through the structure. The accelerating fields in cylindrical coordinates are given by:

c E B k k r k E E k k k r k E E r t z t r z t t z = θ = − = θ = ϕ ) ( I ) cos( ) ( I ) sin( 1 2 0 2 0 [2-9]

The amplitude and phase of the accelerating field due to the rf input power P and the rf phase ϕ are given by:

ϕ + − ω = θ α =

ò

z z s dz z k t P R E 0 0 ' ) ' ( 2 [2-10]

The phase velocity of the accelerating wave increases along the structure in such a way that an electron entering the linac with initial energy _γ0mc2 is accelerated but remains at

the same position relative to the wave. This requires the linac to be operated at its design power PO and the particle must be input at the design phase θ0 with respect to

the ‘crest’ of the wave. At these settings, the longitudinal electric field of the accelerating wave at the position of the design particle is given by:

) cos(

2 0 0 0

0 α R P θ

E = s [2-11]

To calculate the integrated kz in [2-10], we use the particle’s Lorentz factor

(31)

the longitudinal position of the particle relative to the beginning of the structure. This yields:

(

)

2 0 2 0 1 1 ) ( ) ( Fz z z F z + − = Þ + = γ β γ γ [2-12]

where β=v/c and F the normalized design acceleration and:

F E q m c e e = − 0 2 [2-13]

Combining [2-12] with kz=k0/β(z) and k0=ω/c yields:

[

]

k z dz k F Fz z z ( ' ) '= ( + ) − − −

ò

0 0 2 0 2 0 1 1 γ γ [2-14]

The beamloading wave has the same structure as the accelerating wave, i.e. the same wavenumbers, but different amplitude and phase. When the linac is run at its design values, the beamloading field opposes the accelerating field and its amplitude increases along the structure. When the linac is used far off its nominal settings, general analytical expressions or correct approximations for the beamloading wave can not be derived. To allow this element to be used for all input settings, differential equations for the beamloading amplitude and phase are solved while tracing the particles.

The accelerating rf-wave has amplitude E and phase θ as defined in [2-10]. The beamloading wave has amplitude Eb and phase θ+θb. Because the structure of both

waves is similar, the total fields are:

(

)

(

)

c E B k k r k E E E r k E E E r t z t b b r t b b z = θ + θ + θ = θ + θ + θ = ϕ ) ( I ) cos( ) cos( ) ( I ) sin( ) sin( 1 0 [2-15]

To produce a constant accelerating gradient, the power flowing through the constant gradient accelerator decreases linearly as:

) 2 1 ( 2 0P 0 P P 0 0z dz dP z z Þ = − α α − = = = [2-16]

(32)

Because this expression both holds for the accelerating and the beamloading wave, it can be substituted into [2-10] to obtain the relation between boamloading power Pb

power and electric field Eb as function of longitudinal position:

z P R E s b b 0 0 2 1 2 α α − = [2-17]

For a repetitive train of bunches, the steady-state power transferred from the particle beam into the beamloading wave is given by the product of the velocity of the particles in the direction of the accelerating field times the beam current per particle Ib/N:

E v⋅ = N I dt dP b [2-18]

Differentiating [2-17] and combining it with [2-18] directly yields the differential equation for the amplitude of the beamloading wave:

N I z v R dt dEb s z b 0 0 2 1₋ _α α = [2-19]

To calculate the phase of the beamloading wave, it is characterized by the amplitudes u and v of two independent waves, 90° out of phase.

î í ì − = + + = + Þ î í ì = = ) sin( ) cos( ) cos( ) cos( ) sin( ) sin( ) sin( ) cos( θ θ θ θ θ θ θ θ θ θ v u E v u E E v E u b b b b b b b b [2-20]

The differential equations for u and v to be solved are then given by:

) cos( 2 1 ) sin( 2 1 0 0 0 0 θ α − α = θ α − α = N I z v R dt dv N I z v R dt du b z s b z s [2-21]

with boundary condition ut=0=vt=0=0

The u and v representation is converted to power and phase to be written in the GPT outputfile using: ) / arctan( 2 2 u v v u E b b = + = θ [2-22]

The typical evolution of the amplitude and the phase of the beamloading wave in the first FELIX linac is shown in Figure 2-3. The amplitude increases along the entire

(33)

length of the structure and the phase immediately becomes π, which means the wave is precisely in anti-phase with the accelerating wave, as expected.

1 2 3 4 5 z [m] 0 2 4 6 Eb [M V /m ] 1 2 3 4 5 z [m] 0 1 2 3 4 p has e [ ra d ]

Figure 2-3: Amplitude and phase of the beamloading wave for the first FELIX linac when a 0.22 nC bunch with average current of 0.2 A is accelerated from 4 to 25 MeV. The entrance of the linac is at 1.58 m and the exit at at 3.98 m.

2.4.2 Line segment with current

The linecurrent element models a straight line segment, running from (x1,y1,z1) to

(x2,y2,z2), carrying a current I. A number of linecurrents are typically chained together

to model the magnetostatic fields of an irregularly shaped solenoid, see section 5.4.3. Technically, the only two parameters the linecurrent element needs are the total length

L and the current I because the ECS specification allows the element to be positioned

anywhere in 3D space.

To simplify the equations and to save CPU time, a new coordinate system is derived in which the z-direction is parallel to the line segment and the origin is chosen to be the center point of the line:

2 / ) , , ( 1 2 1 2 1 2 line= x +x y +y z +z o [2-23]

The z-column of the coordinate transform matrix of Mline, the projection of the vector

(0,0,1), can be calculated by the normalized direction of the difference between start and end position:

L z z M L y y M L x x M / ) ( / ) ( / ) ( 1 2 33 1 2 23 1 2 13 − = − = − = [2-24]

(34)

The x- and y-columns are not be uniquely defined because of the cylindrical symmetry of the line. However, they need to be both orthonormal to each other and to the z-column. Filling in zero’s for the x and y-columns and applying Singular Value Decomposition is an option, but the following simple recipe is used because of the low dimensionality of the problem:

• The new x-column is set to have a 1 only at the position of the smallest absolute value in the z-column.

• The new y-column is calculated by: y₌(z_×x)/|z_×x| • The x-column is overwritten using: x₌(y_×z)/|y_×z|

The total coordinate system transformation is the product analogous to equation [2-8] of the element coordinate system (M,o)ECS with the line position and orientation matrix

(M,o)line.

In the new coordinate system, the magnetostatic field of the element is given by:

2 1 2 2 1 2 2 1 2 2 / 2 / 2 2 2 0 ) ( ) ( ) ( 1 , 0 , 0 z z y y x x L dz z y x L z L z − + − + − = ÷ ÷ ø ö ç ç è æ + + × ∇ = +

_ò

− π µ 4 I B _[2-25]

where the integral and curl are calculated analytically, resulting in the following equations in cylindrical coordinates:

(

)

(

)

(

)

(

)

÷÷ ÷ ø ö + + + + + + ç ç è æ − − − − + − − = + = ) ( ) ( 1 ) ( ) ( 1 4 2 1 2 2 1 2 2 2 1 2 2 1 2 2 1 2 2 1 2 0 2 2 2 z L z L r z L r z L z L r z L r B y x r π µ ϕ I [2-26]

(35)

2.4.3 Bar magnet

The barmagnet element models a bar shaped magnet with permanent uniform magnetization M. It has been used to model the FEM undulators [40] for example. The magnet is centered around the origin with x-, y- and z-dimensions a, b and L respectively as shown in Figure 2-4.

Figure 2-4: Barmagnet geometry.

The H field produced by the barmagnet element can be derived from the fields of two identical plates located at x <1

2a, y <21b, z= ±12L with a uniform magnetic surface

charge density1:

0

µ

σm=±M [2-27]

The field of a magnetically charged plate at z= 0, x _<1

2a, y <12b is:

(

) (

)

÷÷_ø ö ç ç è æ + − + − − =

_{ò ò}

− − 2 2 2 2 2 2 2 plate ' ' ' ' 1 4 grad a a b b dx dy z y y x x π σ H [2-28]

where the integrals and gradient can be evaluated analytically:

(

)

b y y a x x b y y a x x z r y x x r y r z y x 2 1 2 1 2 1 2 1 , , plate / arctan ) log( ) log( 4 ) , , ( + = + = − = − = ú ú ú û ù ê ê ê ë é ÷ ÷ ÷ ø ö ç ç ç è æ + − + − = π σ H [2-29]

1_{Although we did not test this, an alternative approach would be to model the magnet by a uniform}

(36)

The total field is the sum of the two plates: ) , , ( ) , , ( ) , , ( plate 12 plate 21 bar x y z =H x y z+ L −H x y z− L H [2-30]

The magnetic B field is then given by: magnet the inside ˆ magnet the outside 0 0 z H B H B M + = = µ µ _[2-31]

2.4.4 Solenoid with rectangular cross section

The rectcoil element models a solenoid with rectangular cross section and uniform current density, as shown in Figure 2-5. The field off axis is approximated using a 4th order power series expansion of the analytically calculated field on axis. This basic element has been used for a variety of applications.

Figure 2-5: Geometry of the rectcoil model with inner radius r1, outer radius r2, total length L and current I.

To calculate the field on axis, the magnetic field of a single current loop is integrated in both the z and r direction:

dr dz z r r I r z B a z a z z

ò ò

+ − + = r2 r1 2 3 ) ( 2 ) , ( ₂ ₂ 2 0 µ [2-32]

where a=L /2. The analytic solution to this double integral is:

(

)

) , ( ) , ( ) , ( ) , ( ) , ( log ) ( 4 ) , ( 1 2 1 2 2 2 1 2 0 r a z IBz r a z IBz r a z IB r a z IB r z B z r r z r r a I r z IB z z z z − + − − + − + = + + − = µ [2-33]

(37)

The off-axis field is approximated using a 4th order power series expansion [12]. As a result, the fields are only correct near the axis of the element.

B z r B z B z r B z r B z r B z r B z r z r ( , ) ( ) ( ) ( ) ( , ) ( ) ( ) '' '''' ' ''' = − + = − + 1 4 1 64 1 2 1 16 2 4 3 [2-34]

Calculating the derivatives is tedious when done by hand, but not very difficult. Because the resulting equations are not very illuminating, they are not shown here.

2.4.5 Field maps

Field maps are used in the GPT code to model electromagnetic fields calculated by external field solvers. Various 2D and 3D field map importers are implemented for electrostatic, magnetostatic and TM cavity structures. They are used for a variety of applications, especially where analytical expressions can not be used. All GPT field map elements start by reading a file containing the calculated fields interpolated on a rectangular grid. Various utility programs are available to assist in the conversion from the output of several commercial field solvers like Superfish [13] and Tosca [14]. Alternatively, a tabulated ASCII file can be used to present the data to GPT.

The rows containing the data, for example r, z, Br and Bz for the cylindrically symmetric magnetostatic element, do not need to be in any particular order because the complete set is first sorted for convenient internal use. From this sorted data, the lower and upper bounds and grid spacing are calculated in all directions. Then it is made sure that all the grid points are present and lie within a range of ±10–3 _{∆ off the expected}

position, where ∆ is the grid spacing in the corresponding direction.

Bilinear interpolation is used to calculate the electromagnetic field at a specified particle position. In two dimensions only the nearest four grid points, numbered 1 to 4 as shown in Figure 2-6 are required. The 3D elements require 8 points. Because the grid must be rectangular and all points must be present, the CPU time required to find the nearest points is independent of the size of the field map.

(38)

Figure 2-6: Datapoints around the particle position in 2D (R,z) and 3D (x,y,z).

For a 2D cylindrically symmetric magnetostatic field map, bilinear interpolation yields:

B t u Br t u Br tuBr t uBr B B t u Bz t u Bz tuBz t uBz r z = − − + − + + − = = − − + − + + − ( )( ) ( ) ( ) ( )( ) ( ) ( ) 1 1 1 1 0 1 1 1 1 1 2 3 4 1 2 3 4 ϕ [2-35] where t=

₍

r−r1

₎

∆r u=

₍

z−z1

₎

∆z

A 2D cylindrically symmetric cavity in TM mode requires the angular frequency ω and the rf phase offset ϕ as additional parameters. The fields themselves, again bilinearly interpolated, are then given by:

(

)

(

)

(

1 2 3 4

)

4 3 2 1 4 3 2 1 ) 1 ( ) 1 ( ) 1 )( 1 ( ) sin( ) 1 ( ) 1 ( ) 1 )( 1 ( ) cos( ) 1 ( ) 1 ( ) 1 )( 1 ( ) cos( ϕ ϕ ϕ ϕ ϕ ω ϕ ω ϕ ω ϕ t t u B t u B tuB t uB B uEz t tuEz Ez u t Ez u t t E uEr t tuEr Er u t Er u t t E z r − + + − + − − + − = − + + − + − − + = − + + − + − − + = [2-36]

As a last example, 3D bilinear interpolation for an electrostatic field map is given by:

E t u v Ex t u v Ex tu v Ex t u v Ex

t u vEx t u vEx tuvEx t uvEx

E t u v Ey t u v Ey tu v Ey t u v Ey t u vEy x y = − − − + − − + − + − − + − − + − + + − = − − − + − − + − + − − + − − ( )( )( ) ( )( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( )( )( ) ( )( ) ( ) ( ) ( ) ( )( ) 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + − + + − = − − − + − − + − + − − + − − + − + + −

t u vEy tuvEy t uvEy

E t u v Ez t u v Ez tu v Ez t u v Ez

t u vEz t u vEz tuvEz t uvEz

z ( ) ( ) ( )( )( ) ( )( ) ( ) ( ) ( ) ( )( ) ( ) ( ) [2-37] where t=

₍

x−x1

₎

∆x u=

₍

y−y1

₎

∆y v=

₍

z−z1

₎

∆z

More sophisticated interpolation schemes, for example based on splines, provide higher order smoother fields with significantly less grid points. However, because computer memory is generally not a concern for 2D fields, the more robust and faster bilinear interpolation is used. Because the memory argument does not hold for the 3D case, direct interpolation from the finite element mesh used in the external calculations would under most circumstances be preferable to interpolation from a rectangular grid.

(39)

2.4.6 HE waveguide mode with particle wave interaction

The HEbm element models a single HEmn mode in a corrugated rectangular waveguide

including particle-wave interaction. It is used in section 5.2 and an example of how GPT can be used as an FEL code by solving additional differential equations while tracing the particles. Equations for other waveguide modes can easily be derived following the recipe presented here. The parameters of the HEbm element are listed in Table 2-B.

Table 2-B: Parameters of the HEbm element

parameter Description

ECS Element Coordinate System

a Total wave guide length in x direction [m].

b Total wave guide length in y direction [m].

L Total wave guide length in z direction [m].

P Power flowing in the specified HE mode [W].

ω Angular frequency [s–1].

ϕ Phase factor.

m Mode number in x direction.

n Mode number in y direction.

Lb Bunch length [m].

The electromagnetic fields of a HEmn eigenmode in a corrugated waveguide are given

by: ÷ ÷ ÷ ÷ ÷ ø ö ç ç ç ç ç è æ + − − + − = ) cos( ) sin( ) cos( 0 ) sin( ) sin( ) sin( ϕ ω ϕ ω z k t y k x k k k z k t y k x k A z y x z x z y x E [2-38] ÷ ÷ ÷ ÷ ÷ ÷ ÷ ø ö ç ç ç ç ç ç ç è æ + − − + − + + − = ) cos( ) cos( ) sin( ) sin( ) sin( ) sin( ) sin( ) cos( ) cos( 2 2 ϕ ω ϕ ω ϕ ω ω z k t y k x k k z k t y k x k k k k z k t y k x k k k k A z y x y z y x z z x z y x z y x B [2-39]

where A is a linear multiplication factor for all field amplitudes. The wavenumbers kx, ky and kz are calculated from the mode numbers m, n, the waveguide dimensions a, b

(40)

2 2 2 2 y x z y x k k c k b n k a mπ k − − = = = ω π [2-40]

The total power P flowing in the HEmn mode can be expressed as a function of the field

amplitude A using: z z x g kk k k c ab A v ab w P 2 2 0 2 8 + = = µ [2-41]

where vg is the group velocity vg =ckz/k and the average field energy density is obtained by calculating the average electromagnetic energy in a box with the dimensions of one oscillation of the field:

( )

2 2 2 2 0 0 0 0 0 3 0 8 2 2 2 2 z z x z y x k k k A dV k k k w kx ky kz ₊ = ⋅ = ⋅ =

_{ò ò ò}

ε ε π ε π π π E E E E [2-42]

The interaction between the particle beam and the HEmn wave is derived from

conservation of energy: All energy lost by the electrons due to movement in the direction of the electric field must be gained by the light level.

The power lost by all charged particles qi, due to a decrease in kinetic energy caused by

the electric field of a mode, is given by:

å

⋅ = i i i q dt dW E v lost [2-43]

This lost power must be matched by an increase in radiation power:

2 2 2 2 0 2 light light lost 8 0 z z x b b k k k c L ab A L ab w W dt W d dt W d + = = = + µ [2-44]

where Wlight is the energy of the radiation field in the total bunch length Lb. The

assumptions are made that A does not increase significantly on scales of the order of the bunch length and that the bunch length can be regarded a constant.

The general particle tracer code : design, implementation and application