Accurate modelling of charged particle beams in linear accelerators

A higher order, energy conserving discretization technique for beam dynamics simulations is presented. The method is based on the discontinuous Galerkin (DG) formulation. It utilizes locally refined, non-conforming grids which are designed for high spatial resolution along the path of charged particle beams. Apart from this formulation, the paper introduces a class of general symplectic integrators which conserve discrete energy in a modified sense. Specialized split-operator methods with optimum dispersion properties in the direction of particle motion are, additionally, derived. The application examples given in the paper are performed in a high performance computing environment. They include the self-consistent simulation of the RF electron gun developed by the Photo Injector Test Facility at DESY Zeuthen (PITZ) project and the computation of short range wake fields for ultra-relativistic electron bunches.


Introduction
The development of the numerical methods presented in this study is motivated by the need for accurate beam dynamics simulations in linear accelerators. In particular, the low emittance electron source for Free Electron Lasers (FEL) developed at the Photo Injector Test Facility at DESY Zeuthen (PITZ) [1] serves as a major application throughout the paper. The most critical issue in the simulation of such structures is the extreme multiscale character of the problem. In the case of the PITZ injector, e.g., electron bunches of a few mm length are considered. At the same time, the effects of accelerator geometry over a distance of several meters must be taken into account.
Beam dynamics simulations are typically based on the Particle-In-Cell (PIC) method [2], which involves the coupled solution of Maxwell equations and the particles equations of motion. A widely used technique for the solution of Maxwell equations is the finite integration technique (FIT) [3]. It is often used in conjunction with regular Cartesian grids, although some work has been done to generalize its application to boundary fitted and unstructured grids [4,5]. The most appealing property of FIT is, its mimetic form to Maxwell equations in continuum. In addition, the discrete differential operators used in FIT guarantee by construction strict energy and charge conservation [6]. The former is, in particular for long time beam dynamics simulations, a crucial property, since the conservation of electromagnetic energy is pre-requisite for phase space conservation, and thus, for an accurate beam characterization.
One difficulty in grid-based beam dynamics simulations is the presence of short range space charge fields in high current beams. Capturing the space charge effects in the simulation requires 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT an extremely high spatial resolution for the particle bunches over a large computational domain. In order to achieve such resolutions with affordable computational resources, the discretization grids should be adapted to the problem at hand by providing local mesh refinement. The approach we use in this study is that of non-conforming, octree type Cartesian grid discretization. There are three main reasons for employing these type of grids. Firstly, the high aspect ratio between bunch dimensions and surrounding geometry are well matched with a minimum of discretization cells. Secondly, for PIC simulations, highly efficient search and interpolation algorithms taking advantage of the Cartesian cell geometry are available. Thirdly, the simple topology description of these grids allows for a simple implementation of the numerical code.
The use of non-conforming grids is the motivation for considering a discretization technique based on the discontinuous Galerkin (DG) finite element method (FEM). As in the standard FEM, the method uses a set of basis functions to project Maxwell equations into a finite elements space. However, FEM typically uses globally continuous approximation functions on conforming grids. In contrast, the DG-FEM method employs a discontinuous approximation space. The resulting discrete equations are inherently explicit in time, which makes the method particularly efficient for time domain simulations.
The basic DG-FEM scheme has been thoroughly investigated by Cockburn and Shu [7] in the context of the compressible Navier-Stokes equations. Jacobs and Hesthaven [8] implement a higher order dissipative DG-FEM method for Maxwell equations with a Runge-Kutta integrator. Recently, Canouet et al [9] have proposed an energy conserving DG scheme on non-conforming Cartesian grids. This is also the approach we have adopted in this paper. However, our formulation is based on discrete differential operators, which are derived from the weak form of Maxwell equations. Analogously to FIT, the energy conservation of the scheme can be shown to hold by virtue of the properties of these operators alone. This formulation allows us also to systematically construct higher order symplectic time integrators, so that the approximation orders of both, time and space discretizations, match.
Another critical issue in beam dynamics simulations is numerical dispersion. Due to Lorentz contraction, the electromagnetic field of an ultra-relativistic electron bunch is highly anisotropic. The short wave lengths of the electromagnetic field spectrum concentrate in the direction of particle motion, thus, producing steep field gradients in this direction. This part of the electromagnetic field spectrum, being close to the cut-off wavelength of the discretization grid, propagates at the wrong numerical phase velocity. Apart from employing higher order numerical approximations with possibly local grid refinement, it is possible to construct numerical schemes which specifically address the issue of field anisotropy in beam dynamics simulations. In this study, we present a class of low dispersion split-operator methods which can be used within the frameworks of both, FIT and DG discretizations.
The paper is organized as follows. For the purpose of later reference, we provide a brief review of the PITZ project and of FIT in sections 2 and 3, respectively. In section 4, we introduce the DG formulation for the Maxwell equations and deduce the properties of DG discrete operators leading to energy conservation. Then, the higher order symplectic integrators and the respective modified energy expressions are derived. Also in this section, a general domain decomposition technique for parallel computations with locally refined grids is described. In section 5, low dispersion split-operator methods for use with either DG or FIT discretizations are proposed. Finally, in section 6, two large scale applications of the methods are given. The first is the self-consistent simulation of the RF gun of the PITZ injector. The second, is the computation  of short range wake fields in the ultra-relativistic regime using a longitudinally dispersion free split-operator method and a moving window approach.

The PITZ project
The purpose of PITZ is the development and testing of a high quality, pulsed electron source for use in the European X-Ray Laser Project (XFEL) [10] and in the future International Linear Collider (ILC). Besides the high quality concerning beam emittance, the injected electron bunches have to be highly charged in order to achieve the phase-space density required for the Self-Amplified Spontaneous-Emission (SASE) to occur in the downstream FEL. The extraction of a bunch charge of 1 nC from an emission area of ∼3 mm 2 and a duration of ∼20 ps leads to a space charge dominated beam with local current densities above 1 kA cm −2 . This exceeds the limit which can be provided by thermal cathodes. Therefore, laser-driven photocathodes are used which, additionally, allow for the length and shape of the bunch to be controlled on a ps timescale. Figure 1 shows the geometry and the main design parameters of the PITZ RF gun. The emitted bunch is accelerated in a 1.5-cell L-band cavity providing for an accelerating gradient of 42 MV m −1 at an operation frequency of 1.3 GHz.
A low transverse emittance beam is crucial for FEL operation. In the case of the PITZ injector, a global emittance minimum of ∼1.0 mm mrad has to be achieved. In a second phase of the project, the accelerating gradient of the RF gun will be increased to 60 MV m and a TESLAlike booster will be attached behind the injector at the position of minimum emittance [1]. The expected emittance 'freezing' by acceleration to ultra-relativistic energies, however, has still to be verified experimentally.

FIT
A general framework for the self-consistent simulation of electromagnetic fields and particle dynamics in accelerators is provided by the FIT. FIT-based beam dynamics simulations in the context of the PITZ project have been previously presented, e.g., in [11]. The spatial discretization of FIT utilizes a pair of dual-orthogonal staggered grids (G,G). Denoting by ( e, h) the electromagnetic degrees of freedom and by j the coupling current associated with charged particles, the semidiscrete FIT equations read d dt The discrete fields ( e, h) have the precise interpretation of electric and magnetic voltages along the edges of G andG, respectively. The curl-operator C has the structure where (P x , P y , P z ) denote antisymmetric, discrete derivative operators. Finally, the symmetric positive definite material matrices (M , M µ ) provide constitutive relations, for the fluxes of dielectric displacement, d, and magnetic induction, b, defined on the faces of the grid doublet G andG, respectively. A commonly used time stepping method for the solution of (3.1) is the leapfrog integrator. It leads to the conditionally stable time-discrete update equations where the scalings, C 1 = M −1 C, C 2 = M −1 µ C T and j n = M −1 j n are used. In [6], it is shown that, for a simulation domain enclosed by perfectly conducting boundaries, equation (3.1) conserves electromagnetic energy, In the following, we develop a formally similar discretization method which uses locally refined grids and higher order spatial approximations, while maintaining the mimetic structure of (3.1) and the energy conservation property (3.5).

The weak DG formulation
In order to obtain a higher order scheme for the solution of Maxwell equations, we proceed as follows. Given a discretization of the computational domain with non-overlapping cells, C j , for j = 1, 2, . . . , N, we seek an approximation for the electromagnetic field components in the form E(r, t) ≈Ẽ(r, t) = j,p e j;p (t)ϕ j;p (r), where (Ẽ,H) are the approximated fields, (e j;p , h j;p ) are the electromagnetic degrees of freedom and ϕ j;p , for p = 1, 2, . . . , P, denote a set of linearly independent basis functions in each cell C j . The basis functions are required to be continuous within C j , and vanish otherwise. Thus, a Pth order local approximation is obtained, e.g., for the choice where ϕ p (r) are polynomials of order p. However, the approximations (4.1) will in general be discontinuous at the interface between two neighbouring cells, as opposed to standard FEM, where global continuity of at least the tangential field components is imposed (cf [12]). Following the standard Galerkin procedure, we replace (4.1) into the Maxwell equations and test the resulting equations with each of the basis functions. Omitting, for simplicity, the source term as well as explicit time and space dependencies, this yields the weak DG formulation as . . , N, and ∀ q = 1, 2, . . . , P. Equations (4.3) define a set of ordinary differential equations (ODEs) for the degrees of freedom, (e j;p , h j;p ). Note, that the discontinuity of approximations (4.1) does not allow for a straightforward application of Stokes theorem for transforming the second volume integral in (4.3) into a surface flux integral. Instead, the approach used by the DG scheme involves the so-called numerical fluxes. They are defined by where n is the outward directed normal to ∂C i and (Ẽ * ,H * ) are effective electric and magnetic fields at the cell interface. We chose the numerical fields as where (Ẽ − ,H − ) are electromagnetic fields computed according to (4.1) from the inside of C i and (Ẽ + ,H + ) denote fields computed from the inside of a neighbouring cell sharing an interface with C i . At a perfectly conducting boundary interface the numerical fields are imposed to This choice of numerical fields is not unique (see [8] for an alternative approach). Below, however, it is shown that the centred flux approach (4.5) guarantees the strict energy conservation of the numerical scheme.
where (e, h) denote the vectors of electromagnetic degrees of freedom, and M ε , M µ are tensor mass operators of the form The latter are intrinsically block-diagonal and symmetric positive definite for an arbitrary choice of linearly independent basis functions. Indeed, the matrix elements of, e.g., M x ε are explicitly found to which verifies the above properties when a homogeneous material distribution within each discretization cell is assumed. A closer inspection of the curl-operator C appearing in (4.7), reveals that it can be written in the form where (P x , P y , P z ) are easily interpreted as discrete derivative operators. In terms of matrix elements, they are explicitly given by Each of the above equations contains a pure block-diagonal part, which is responsible for the coupling of the cell-internal degrees of freedom. The second term contains the coupling to neighbouring cells, which occurs through interface fluxes alone. This local character of the discrete operators in the DG formulation makes the method particularly well suited for parallel implementations (see subsection 4.5).

8
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Energy conservation
The energy conservation of the DG formulation can be proven using the properties of the discrete operators (4.10) and (4.11a)-(4.11c). The basic property of the derivative operators is their antisymmetry. Indeed, performing an integration by parts in (4.11a) and using n x | ∂C i = −n x | ∂C j for every pair of neighbouring cells, C i , C j , sharing a common interface, it follows in an element by element fashion Thus, the following relations hold true The semidiscrete equations (4.7) can, then, be written as d dt which is a system of ODEs with a pure antisymmetric coefficient matrix. Note that (4.14) are formally the same as the semidiscrete equations (3.1) obtained in FIT. Although, the unknown vectors of electric and magnetic fields have a different interpretation than in FIT and the discrete operators are defined in the higher order discontinuous finite element space (4.2), this analogy proves useful in deriving the energy conservation property of the scheme.
Using the positive definiteness of the mass operators (M ε , M µ ), their square root can be taken, so that (4.14) transforms to d dt where (ẽ,h) andC are defined bỹ In full analogy to FIT [6], the semidiscrete energy is then defined as (4.17) which is clearly the discrete analogon to the electromagnetic energy in continuum. Furthermore, since the coefficient matrix in (4.15) is antisymmetric, The energy conservation of the method means, in particular, that equation (4.14) has a stable solution. Thus, an appropriate ODE integrator can be found, such that stable time stepping upon the time discrete equations can be performed.

9
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Higher order symplectic time integration
The energy conservation property (4.18) reveals that the semidiscrete equation (4.14) preserves the Hamiltonian structure of the Maxwell equations in continuum. For the numerical integration of such systems, the so-called symplectic integrators are widely used. The reader is referred to [13] for a very complete discussion of these integrators in the context of finite dimensional Hamiltonian systems. Recently, higher order symplectic integrators have been used also in discrete electromagnetics for a specific vector finite element type discretization of Maxwell equations [14]. However, the benefits of applying them with continuous FEM are severely limited by the presence of non-local mass matrices. In this case, either the solution of a sparse system of equations in every time step, or a mass lumping procedure is required. The former strongly affects the efficiency of time domain simulations and the latter will generally deteriorate the order of approximation of both, space and time discretizations. We derive the higher order symplectic integrators using an operator splitting approach. For this purpose, a compact notation is introduced by rewriting (4.14) as The formal solution of (4.19) within a single time step t is given by where G( t) is the (unitary) propagation matrix, y n and y n+1 denote the discrete field solutions at the present and next time levels, respectively. The lowest order symplectic integrator is obtained by applying the splitting which is referred to as Godunov splitting [15]. Selecting the split-operators as the upper and lower parts of H, respectively, and using the observation that (H L , H U ) are nilpotent operators, the splitting (4.21) leads to the approximation Then, in terms of discrete field variables the first order time update equation is (4.24) Apart for a time shift of t/2 in the initial value ofẽ, the above scheme is formally the same as the leapfrog integrator (3.4) used with FIT. The curl-matrices appearing on the righthand side of (4.24), however, are specific DG operators corresponding to higher order spatial discretizations. Also, (4.24) represents purely explicit time stepping, regardless of the order of basis functions, since the block-diagonal mass matrices M ε , M µ can be trivially inverted (see also subsection 4.4). A second order time update scheme is obtained by applying the Strang splitting procedure [16] on the formal solution (4.20): so that (4.26) The full update equation in this case is given by It can be easily verified that (4.27) is identical with the Verlet method, which is commonly used in the solution of N-body problems. Third and fourth order update schemes can be constructed more systematically by employing the Suzuki decomposition formula for exponential operators [17], so that, using the splitting (4.22), the propagation matrix is approximated to In particular, the choice yields the third order Ruth's scheme [18] and gives the fourth order scheme of Candy and Rozmus [19]. Both schemes can be implemented by sequentially applying the first order update (4.24) with the respective scaling coefficients (α i , β i ). Thus, no additional storage is required for higher order time integration.
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Equations (4.24) and (4.27)-(4.31) provide a complete platform for obtaining up to fourth order accurate discrete solutions of Maxwell equations in both, the space and time variables. We will postpone the discussion on the stability of the proposed schemes to a forthcoming paper. We only note, that for each of the update equations employing DG spatial operators, conditional stability can be proven. In particular, in the simplest case of discretization with piecewise constant basis functions and for a regular Cartesian grid of spacings ( x, y, z), it has been shown in [20], that (4.24) is stable under the Courant condition, c t 2/ x −2 + y −2 + z −2 . Considering the discrete electromagnetic field energy (4.17), it can be verified, e.g., using the update (4.24) that for two subsequent time steps The energy of the time discrete fields is, thus, no longer conserved. However, according to [13], it is possible to define a modified Hamiltonian for the symplectic mappings of the form (4.28) corresponding to a time discrete conserved energy. Instead of deriving the modified Hamiltonian from the electromagnetic field update equations, here, we follow a simpler approach based on algebraic arguments for defining the modified energy of schemes (4.24) and (4.27).
In the case of the first order scheme and using the notation of (4.19), we rewrite the update equations as remains constant in every time step. It is a first order correction to the semidiscrete energy (4.17) and, if the time step is small enough for W (1) to be positive, it can be interpreted as a time discrete electromagnetic energy. This requirement provides at the same time a sufficient condition on t for the stability of (4.24). Similarly, the second order scheme (4.27) can be transformed to In this case, the modified discrete energy is given by TCTCh (4.37) and, if W (2) > 0, it is a second order correction to the semidiscrete energy (4.17).
In figure 2, the overall accuracy of the method is tested in the simple case of a EM 111 mode in a box resonator with uniform grid discretization. The numerical error for different orders of  spatial approximation using the first order update (4.24) is shown in figure 2 (left panel). The dispersion error of the numerical scheme decreases drastically with increasing approximation order. Figure 2 (right panel) shows the behaviour of discrete energy, again, for different orders of spatial discretizations. The discrete energy shows no secular increase or decrease with time. However, it is not strictly conserved. In contrast, the modified discrete energy is verified to remain constant.

Local mesh refinement with non-conforming grids
Non-conforming grid discretizations are constructed by recursively subdividing the cells of an initially regular Cartesian grid into smaller ones in the regions which are expected to be critical for the electromagnetic field solution. These include, e.g., the path of the electron bunch as well as the metallic edges of accelerating cavities. The procedure results in an octree data structure containing cell geometry, refinement level, parallel decomposition information (see subsection 4.5) and the field data storage associated with each of the cells. The grid connectivity is explicitly stored in order to efficiently use the discretization data for particle localization and interpolation. Figure 3 shows a transversal slice of such a discretization for the PITZ gun geometry. Although the grid has been produced for illustration purposes only with less than 250000 active cells, it already provides for a maximum resolution of ∼250 µm. This corresponds to an aspect ratio of 1/256 between the largest cells edge outside of the gun geometry and the finest discretization level along the path of the electron bunch. For the non-conforming Cartesian grids presented here, the formulation (4.3) simplifies considerably since a tensor product basis of orthogonal polynomials can be used. The basis functions used in our implementation are,

38)
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 3. Example discretization of the PITZ RF gun using a non-conforming grid.
where i = 1, 2, . . . , N denotes the cell index, (p, q, r) = 1, 2, . . . , P are the spectral indices of the 1D-approximation basis and L p (x) is the Legendre polynomial of order p. For this choice, the mass operators (4.9) reduce to purely diagonal material matrices.

A flexible domain decomposition method for parallel computations
The parallelization model used in the simulations is based on the distribution of computational tasks and simulation data among a number of memory-independent processes in a computing cluster. Simulation efficiency, however, may only be achieved if interprocess communication is kept at the lowest level possible. Additionally, a well-balanced workload should be assigned to each of the processes in order to ensure coherency in parallel program execution. For these purposes, we use a geometric decomposition of the computational domain into subdomains which are, then, assigned to the single processes. In this approach, each process is responsible for (i) performing the electromagnetic time update equations for discretization cells contained in the subdomain, (ii) integrating the equations of motion for numerical particles within the bounds of the subdomain and (iii) exchanging particle and field data with processes assigned to the neighbouring subdomains. The domain decomposition approach is shown schematically for a three-node cluster in figure 4(a). Starting with the whole computational domain, we apply an orthogonal recursive bisection of the domain bounds. This procedure results in a binary tree structure, whose internal nodes are intermediate subdomains whereas the leaf nodes correspond to the active (computational) subdomains. Because of the local nature of discrete operators (4.11a)-(4.11c), only field degrees of freedom residing at active subdomain boundaries need to be exchanged between the processes. The orthogonal decomposition approach minimizes the number of such boundary cells and, therefore, the communication overhead in the field update equations.
In order to determine the optimum decomposition which results in equally balanced workloads among the processes, we follow the scheme proposed in [21] for time domain FIT simulations with structured Cartesian grids. The recursive bisection procedure is performed on  the basis of computational weights w i which are assigned to each discretization cell C i . The total computational load associated with an intermediate subdomain is, w = w i , where the summation includes only cells within the subdomain. If the subdomain has to be decomposed between N processes, the bisection bounds are chosen such that where (w left , w right ) and (N left , N right ) are the computational weights and the number of processes, respectively, assigned to the two subdomains created by subdivision. The algorithm allows for an almost ideally balanced distribution of computational workloads, even for highly refined grids. In addition, it can be applied to simulations involving an arbitrary number of processors. Figure 4(b) shows the resulting decomposition of the gun section of the PITZ injector on a cluster of ten processors.

A class of specialized split-operator methods
The basic idea of split-operator methods is to modify the time domain update equations of the type (4.15), such that certain spatial directions are handled with a greater accuracy than the others. In this manner, the numerical dispersion error in the direction of bunch motion is reduced and, under certain conditions, completely eliminated. We note, that as far as these methods are concerned, the curl-operators appearing in (4.15) are arbitrary discrete operators derived by a consistent spatial discretization scheme such as FIT or a higher order DG formulation.

Longitudinal-transversal (LT) splitting
Decomposing the curls into a part containing only longitudinal derivatives, C l , and a part containing only transversal derivatives, C t , equations (4.19) take the form with a formal solution within a time step t given by Following the same procedure as in section 4, a second order accurate approximation in time of (5.2) is obtained by applying the multiplicative Strang splitting It is readily seen that, each of the above exponentials describes the evolution of a system of the form (4.15), where in particular, e −H l t is the time evolution operator of a 1D-plane wave propagating in the longitudinal direction. The discrete time update of the LT scheme is obtained by replacing each of the evolution operators in (5.3) with the corresponding propagation matrix of the Verlet integrator where G t and G l are further decomposed into lower and upper parts according to (4.26).

Plus-minus (PM) splitting
An alternative split-operator scheme is obtained by using the decomposition whereC p andC m contain only positive and negative derivatives, respectively. Then, we apply a second order additive Strang splitting [22] for approximating the formal solution of (5.5) within a time step by

Dispersion properties of the split-operator schemes
Numerical dispersion relations for the split-operator schemes using the spatial operators of FIT have been derived in [23]. Here, we only show the behaviour of the numerical phase velocity of the LT-and PM-schemes as compared to the leapfrog integrator. Figure 5 shows points per wave length) for waves propagating in the xz-plane have been used. Figure 5(a) shows, that the dispersion error of the leapfrog scheme attains its maxima in the directions of coordinate axes. Contrary, the LT-scheme in figure 5(b) shows a minimum dispersion error in the longitudinal z-direction, and the PM-scheme in figure 5(c) minimizes dispersion in the direction of both coordinate axes. Thus, the effect of the above splittings consists in rotating the optimum dispersion directions in the longitudinal direction. Additionally, the split-operator schemes with homogeneous grids are stable for the 1D-Courant time step t = x/c. This yields a maximum stable time step which is √ 3 larger than the one obtained with the standard leapfrog integrator (cf [23]). For this particular choice of the time step, both split-operator schemes have exact propagation along the z-direction. The dispersion free integration of longitudinal waves in the split-operator approach will be used in subsection 6.2 for wake field computations on moving grids.

Self-consistent simulations of the PITZ injector
PIC simulations were performed for the RF gun of the PITZ injector and the bunch characteristics were extracted. In particular, the bunch emittance and RMS-radius are largely shaped in this section and cannot be easily controlled in the latter accelerator sections. Two discretization techniques, the DG formulation with quadratic approximation basis and FIT with LT-splitting were used in the simulations. From experience it is known that, in order to account for the high space charge fields at the gun cathode, a mesh resolution below 20 µm is needed (the RF gun is 18.5 cm long). This resolution was achieved in the DG simulations using the local mesh refinement of the type shown in figure 3.
The simulation results for a 1 nC bunch produced by a laser pulse of duration 22 ps and RMS-radius of 0.5 mm are shown in figures 6 and 7. In figure 6 (left panel) the RMS-bunch energy for the first 35 cm of the section is shown. Figure 7 shows the transversal bunch emittance and RMS-radius over the same distance. In all cases, the results obtained using the two types of discretization agree extremely well. Due to the local mesh refinement, only about 5.4 million discretization cells were used in the DG simulations as compared to ∼80 million needed by the structured-grid FIT simulations. However, the number of unknowns in the DG formulation is higher, so that the overall performance efficiency for both types of simulations was comparable.

Computation of short range wake fields
Behind the injector section, the electron bunch is accelerated to nearly the speed of light in vacuum. Therefore, electromagnetic fields scattered by a geometrical discontinuity of the accelerator structure will catch-up with the bunch at a very long distance (time) behind the discontinuity. In these situations, integral measures such as the longitudinal wake potential are of interest. In (6.1), Q is the bunch charge and s denotes the longitudinal distance relative to the bunch centre. The computation of wake fields for long accelerator structures is generally accessible only to numerical simulations. In [24]- [26], 2D-simulation codes for rotational symmetric geometries have been presented. Recently, Zagorodnov and Weiland [27] introduced a semi-implicit method for 3D-time domain wake field simulations. In this study, we show simulation results obtained with the code Parallel Beam Cavity Interaction (PBCI) which is designed for massively parallel wake field computations in arbitrary 3D-geometry. The algorithms used include the domain decomposition approach of subsection 4.5 and the purely explicit LT split-operator scheme of subsection 5.1 with either DG or FIT spatial discretizations. In particular, the dispersion free integration in the longitudinal direction allows for a moving window implementation [28]. The latter enables long distance wake field computations with a minimum on computational effort, since the discretization window encloses the bunch region only. Note, that this approach cannot be used with standard time domain methods, because the numerical phase velocity of longitudinal waves must exactly match the speed of light in vacuum. Figures 8 and 9 show the simulation results for a 9-cell superconducting accelerator cavity of the TESLA XFEL project at DESY, Hamburg. The wake potentials and electromagnetic fields of a Gaussian bunch of electrons with a total charge of 1 nC and a RMS-length of 5 mm. The simulation was performed over a total distance of 1.5 m using the moving window technique and a DG spatial discretization with piecewise constant basis functions. More than 80 million discretization cells on a regular Cartesian grid were needed for accurately resolving geometry and bunch extension.
The second application considers the wake field effects of full 3D-geometry in the diagnostics cross of the PITZ injector. In figure 10 (right panel), the geometrical layout of this section is shown. Our investigation consisted of three separate simulations for performing a comparative study of the wake fields produced by the different geometrical components of the structure. In the first simulation, the geometry was simplified to the beam tube including the small kink at the entrance of the section. The second simulation included the cross resonator without shielding tube. The third simulation considered the full geometry as shown in figure 10. The simulation results for an electron bunch of charge 1 nC and RMS-length of 2.5 mm are shown in figures 10 and 11. For resolving the small geometrical details of the structure, a total of 250 million cells distributed on a 24-node cluster of computers were used in the discretization. It was found, in particular, that the kink of 1 mm height within the beam tube is responsible for 10-15% of the wake fields monitored at a distance of 2 m behind the diagnostics cross. Figure 11 demonstrates how wake fields are generated by the geometrical obstacles in this section.

Conclusions
In this paper, we present a collection of novel numerical methods for beam dynamics simulations. We show that the high spatial resolution required by these simulations can be achieved by using locally refined, non-conforming Cartesian grids. The associated discretization technique, the DG-FEM, supports this type of grids and, more importantly, it conserves electromagnetic energy. This property was shown using a discrete operator representation and exploiting the analogy to FIT. Furthermore, higher order symplectic integrators were introduced, which maintain the energy conservation property for a modified, time discrete energy. We derived explicit expressions for the modified energy of the first and second order time updates. We skipped the discussion on the stability and dispersion properties of the method. Instead, issues related to the concrete applications, such as the domain decomposition algorithm for parallel computations and the low dispersion split-operator schemes, were addressed. In particular, the latter enabled us to perform wake field computations in 3D-geometry with the moving window approach.