Extension of the Finite Integration Technique including dynamic mesh refinement and its application to self-consistent beam dynamics simulations

An extension of the framework of the Finite Integration Technique (FIT) including dynamic and adaptive mesh refinement is presented. After recalling the standard formulation of the FIT, the proposed mesh adaptation procedure is described. Besides the linear interpolation approach, a novel interpolation technique based on specialized spline functions for approximating the discrete electromagnetic field solution during mesh adaptation is introduced. The standard FIT on a fixed mesh and the new adaptive approach are applied to a simulation test case with known analytical solution. The numerical accuracy of the two methods are shown to be comparable. The dynamic mesh approach is, however, much more efficient. This is also demonstrated for the full scale modeling of the complete RF gun at the Photo Injector Test Facility DESY Zeuthen (PITZ) on a single computer. Results of a detailed design study addressing the effects of individual components of the gun onto the beam emittance using a fully self-consistent approach are presented.


I. INTRODUCTION
Memory consumption and CPU time represent the main limitations for large-scale electromagnetic field computations. This is especially the case for accelerator physics simulations involving self-consistent charged particle models based on the so-called Particle-In-Cell (PIC) method [1]. Simulations of this type are an indispensable tool for the design and optimization of particle accelerators since they offer a full insight into the beam dynamics down to the particle level. This is especially important for the simulation of low-energy sections of an accelerator, where space charge forces heavily influence the beam but also for, e.g., dark current or electron cloud simulations. Apart from the solution of Maxwell's equations on a discrete grid space, PIC simulations also include the self-consistent solution of the equations of motion. Thus, they provide a full description of the accelerator structure including space-charge and wakefield effects [2].
In this article we address the important problem of numerical efficiency of such beam dynamics simulations using the Finite Integration Technique (FIT). The FIT has been successfully applied for the simulation of a wide range of electromagnetic problems in accelerator physics [3,4]. We propose an extension of this method including dynamic mesh refinement in order to locally adjust the spatial grid resolution according to the dynamics of particles and fields. This leads to considerable savings in the overall number of computational degrees of freedom and, thus, reduces the computational burden in PIC simulations.
The article is organized as follows. After describing the governing equations in Sec. II, a brief review of the FIT on static grids is given in Sec. III, which also serves for introducing the notation. Sec. IV describes the dynamic mesh refinement procedures and Sec. V addresses specific issues of charged particle simulations on dynamically refined meshes. Sec. VI contains two applications. First, an example with known analytical solution is considered for investigating accuracy and efficiency of the proposed method. After, results of a detailed design study of the PITZ photo injector using self-consistent PIC simulations are presented.
The achievements are summarized in Sec. VII.

II. FORMULATION OF THE PHYSICAL PROBLEM
The electromagnetic part of the physical problem considered is described by Maxwell's equations. Their integral form reads where A and V denote arbitrary surfaces and volumes. The magnetic field strength is indicated by H, B denotes the magnetic flux density, E the electric field strength, D the dielectric flux density, and J the electric current density. The electric and magnetic quantities are related according to the constitutive equations For the linear, isotropic, nondispersive materials considered here the permittivity ǫ and the permeability µ are material dependent constants. The conservation of charge follows from Ampère's law (1) and Gauss' law (3) in the form of the continuity equation It reflects that a temporal change of the total charge, Q = V ρ dV , contained in a volume is caused only by a flow of electric current into or out of the volume.
The link between Maxwell's equations and the motion of charged particles is established by the Lorentz force on the one hand and Newton's law on the other hand. Above, q is the charge of the particle, γ is the relativistic factor and u ′ = γ v, with v the velocity of the particle. The mechanical momentum is denoted by p.
The equations of motion, thus, read where u is the position vector and m 0 the particle mass at rest.

III. THE FINITE INTEGRATION FRAMEWORK
For introducing the notation as well as for completeness this section comprises a short review of the Finite Integration Technique. For the spatial discretization of Maxwell's equations, it utilizes a staggered, dual orthogonal grid doublet (G, G) consisting of N p primary nodes, which covers the domain of interest Ω. The staggered grid arrangement is depicted in Fig. 1. Throughout this article, we employ the FIT on Cartesian grids in three-dimensional space. The following discrete integral state variables are introduced where c and A belong to the set of edges and faces of the primary grid (black), and c, A and V belong to the set of edges, faces and volumes of its dual (gray). The quantities ⌢ e and Summing up the voltages ⌢ e and ⌢ h along the four edges enclosing one face yields a discrete but exact representation of Ampère's and Faraday's law, (1) and (2), in the form for every face of the computational grid. The orientation of the involved voltages ⌢ e j , ⌢ h j with respect to the orientation of the loop integral in (1) and (2) determines the summation signs (see Fig. 1).
for every cell of G and G corresponding to the Gauss' law for electric and magnetic charges, respectively.
The constitutive equations (5) and (6) translate to where M ǫ and M µ provide an appropriate mapping of grid voltages into fluxes [3].
The summation signs are gathered together in the matrices C, C and S, S. These matrices are, hence, purely topological. They represent a discrete curl and divergence operator of the FIT. It can be shown that these operators maintain the properties of their counterparts in continuum [5] in the sense that This allows, e.g., for the derivation of a discrete continuity equation [cf. Eqn. (7)] in the

IV. DYNAMIC MESH REFINEMENT
We will restrict the discussion to conformal mesh refinements, which maintain the conformity and duality of the Cartesian grid doublet. A refined mesh is obtained by a sequence of bisections applied to a given mesh cell. The number of consecutive bisections of a cell is referred to as its refinement level L (see Fig. 2).
When mesh refinement is applied the discrete quantities assigned to the primary and the dual grid need to be recomputed according to the modified mesh. In Fig. 3 the refinement of a primary grid cell and the subsequent arrangement of electric grid voltages is illustrated in a two-dimensional cut view. The preservation of the duality of the staggered grids requires adding new nodes to the dual grid, as well as shifting existing nodes. This is shown in Fig. 4. In order to obtain a value for the new or shifted voltages an interpolation procedure has to be carried out. Furthermore, it has to be distinguished between the interpolation of voltages oriented in parallel to the refinement [e.g. ⌢ e ′ x (2 ′ ) in Fig. 3(b)] and those oriented perpendicularly [e.g. ⌢ e ′ z (2) and ⌢ e ′ z (2 ′ ) in Fig. 3(c)]. In the following, two interpolation procedures, linear and high order spline interpolation, will be discussed.

A. Linear Interpolation of Grid Voltages
The application of a linear interpolation for the determination of new or shifted voltages is straightforward. The interpolation procedure described below refers to the refinement scenario shown in Fig. 3 and 4. The lengths of the primary and dual grid edges [cf. (12)] are denoted by |c| and | c |, respectively. Discrete field quantities associated with a specific node i or i are denoted as, e.g., c(i). Primed quantities are either new or they have to be recomputed during the adaptation procedure.  and (c). In addition, the two new nodes 2 ′′ and 6 ′′ have to be inserted. In (b) the assignment of x-oriented magnetic voltages is shown and in (c) the assignment of z-oriented voltages.

a. Refinement
The values of new electric grid voltages oriented parallel to the refine- The splitting of voltages oriented perpendicularly to the refinement [e.g. ⌢ e z (2)] is con- This condition imposes the conservation of the total voltage between the nodes 2 and 3. A detailed sketch of this refinement is shown in Fig. 5. In order to assign values to the refined voltages ⌢ e ′ z (2) and ⌢ e ′ z (2 ′ ) the voltage gradient, i.e., the average electric field, along the axis is approximated by a central difference using the neighboring voltages ⌢ e z (1) and ⌢ e z (3). The voltages assigned to the refined edges read where the sampled electric field variables e z (i) are introduced as the z-coordinate of the node i is denoted by z(i), and the ∆-notation denotes a difference of the respective quantity, e.g., ∆z( 3, 1) = z( 3) − z( 1). Since cells are always divided into halves, |c ′ (2)| equals |c ′ (2 ′ )| and the condition (25) is fulfilled.
The interpolation of the magnetic voltages is more cumbersome. Besides the insertion of new dual grid nodes, also existing nodes have to be shifted in position in order to preserve the duality of the staggered grids. Shifting a node, however, implies a modification of the grid voltages along the associated edges.
The interpolation of magnetic voltages oriented in parallel to the refinement is illustrated in Fig. 4 .
perpendicularly to the refinement, is constrained by the conservation of the total voltage.
However, for a complete description of the refinement algorithm, the bisection of at least two neighboring cells has to be considered. This is illustrated in Fig. 6. For this case, the conservation condition for the total voltage between the nodes 1 and 4 can be expressed as The refined magnetic voltages read with the sampled magnetic field h z ( i) defined as Since the following equations hold true and  is performed The magnetic coarse grid voltages, oriented perpendicularly to the refinement, are obtained by the summation of fractional voltages given on the fine grid (see Fig. 6). The merging of the fine grid magnetic voltages has to obey the constraint (31). The coarse grid voltages are given by Since the following holds true the condition (31) is fulfilled also for the case of grid coarsening.
the coarsening process of the computational grid. In combination with the refinement rules given above, they provide a consistent framework for performing time-adaptive conformal grid refinement within the FIT. Refinement schemes other than by cell bisecting are generally possible. However, these are associated with a larger complexity in mesh administration and will not be considered in this work.

B. Spline Interpolation
A linear interpolation of the discrete field quantities is easy to implement and fast in code execution. The interpolated quantities are second order accurate, which is in agreement with the theoretical order of accuracy of the FIT. However, applying higher order interpolating functions may be beneficial for obtaining a smoother representation of high-frequency fields on the refined mesh.
Polynomials offer a possibility for performing higher order interpolations. However, polynomials of high degrees tend to exhibit an oscillatory behavior and possibly large overshoots, known as Runge's phenomenon [8]. Commonly, spline functions are employed in order to avoid this behavior.
A spline S is a piecewise defined polynomial function of order P which is P − 1 times continuously differentiable, i.e., S ∈ C P −1 [9]. The construction rules determine the specific type of spline. If the spline S is demanded to pass exactly through the given data points, and to be twice continuously differentiable (i.e. S ∈ C 2 ) with the second derivative equal to zero on every interval boundary, the so-called natural cubic spline (C-spline) is obtained.
For its determination, a tridiagonal system of equations has to be solved. During a timedomain simulation adopting time-adaptive grid refinement, hundreds to thousands of grid adaptations are performed, each involving thousands of cells. Thus, the solution of a system of equations is computationally very expensive. In addition, the natural cubic spline may still exhibit overshoots (see Fig. 7).
In order to further mitigate Runge's phenomenon, the conditions on the continuous differentiability of the spline have to be reduced. A spline S of order P which is at most P − 2 times continuously differentiable is called a broken spline or subspline [9]. The Akima spline is a cubic subspline which is an element of the functional space C 1 [10]. Cubic (sub-)splines at every interval boundary, uniquely describes the broken cubic spline function with the coefficients However, the slopes s(i) are, in general, unknown. The Akima procedure makes use of a local heuristic method for the estimation of the slopes. It involves the data point i and two neighboring points on each side. First, the piecewise gradients g, given by are evaluated. These are weighted with the factors w yielding the estimated slope at data point i The fundamental idea of the heuristic can be summarized as: the larger the difference between the two gradients on one side of point i, the larger is the weighting factor applied to the gradient on the other side of this point. This is intended to minimize overshooting and oscillatory behavior. Note that, Akima's heuristic is closely related to the idea of slope limiters. These are commonly used in Finite Volume methods for reconstructing a piecewise continuous solution on the grid (cf. [11]). Using, e.g., the minmod limiter the estimated slope reads Out of its two arguments, the minmod operator chooses the smaller one if they are equally signed and zero otherwise. An overview of slope limiting techniques is found in [12].
In Fig. 7 the interpolation of a set of data points using linear interpolation, the C-spline and the Akima spline is illustrated. The C-spline shows an oscillatory behavior and strong overshooting. For the determination of the Akima spline up to data point 10, the slope was  calculated using Akima's slope definition (50). For the points 11 to 16 the minmod limiter was applied, which effectively avoids any overshooting.
The actual spline interpolation procedure of the discrete field quantities is illustrated in

V. CHARGED PARTICLE SIMULATIONS ON ADAPTIVE GRIDS
Simulations of the dynamics of charged particles using the FIT are performed routinely and have first been reported in 1988 in [15]. For self-consistent simulations the particlein-cell (PIC) algorithm is applied [16]. The PIC algorithm consists of three steps. First, the electromagnetic field at the spacecontinuous position of each macro particle has to be obtained from the space-discrete FIT quantities. This is done by means of a trilinear interpolation. Next, the equations of motion (10), (11) are integrated using the algorithm described in [17]. Finally, the convective currents are calculated from the position increment of all particles. For this last step we adapted the cloud-in-cell (CIC) approach [18,19] to work with time-adaptive mesh refinement. The particles are modeled as a uniformly charged volume of finite extent, i.e., a charged could.
The extent of the cloud is usually chosen to coincide with the size of the cells. However, if a nonequidistant grid is employed or in the case of adaptive grid refinement, the nonuniform sizes of the cells demand for a modification of the algorithm. There are two modification options: first, a constant size of the particle cloud is chosen independently from the grid cell size or, second, the size of the cloud is adapted in order to match the sizes of the involved cells. Depending on the local degree of grid refinement, the first option can largely increase the number of cells affected by one cloud. Besides the coding efforts coming along with this non-local operation, the deposit of fractional currents to a large number of grid points increases the computational load. Therefore, the second option has been pursued.
The adaptation of the cloud in dependence of the grid cell size is illustrated in Fig. 9 for a one-dimensional example.
The modified CIC approach maintains its charge conserving property, meaning that it fulfills the discrete continuity equation (23) in a cell-wise manner.  9. Illustration of the adaptation of the particle cloud size in one dimension. The size of the particle cloud relates to the extent of the covered cells on the dual grid. If the particle is centered within a dual grid cell (cases 1, 3, 5, 7 and 9) the sizes of the cloud and the dual volume coincide.
Otherwise, the cloud has to be adapted asymmetrically around the particle position. For the cases 2 and 4 the charge density within the part of the cloud within the dual volumes 2 or 3 ′ has to be increased as a consequence of the diminished size. Contrarily, for the cases 6 and 8, the charge density of the cloud in the dual volumes 3 ′′ and 4 has to be decreased.

VI. APPLICATIONS
In this section, we present results of the application of the dynamic mesh adaptation algorithm introduced above to the simulation of two setups.

A. Bunch drift in a pipe
We consider a hollow, perfectly conducting pipe, which is closed by a perfectly conducting  Table I. In [20] the analytical solution of this problem for a semi-infinite pipe is given. In the simulations the axis of the cylinder is aligned with the z-coordinate. In the with E z = E z (z(1)), ..E z (z(N z )) T and E z the longitudinal component of the analytical solution of the electric field evaluated at the grid point positions z(i z ).
The total variation is a measure for the smoothness of a function [21]. It is defined as In [22] it is shown that the TV of a function equals the sum of local minima and maxima, where the values at the integration endpoints count once and all other extrema count twice.
The TV is, therefore, a direct measure for oscillations and their amplitudes. For discrete solutions it is computed as The standard FIT is a very well established method for performing PIC simulations.
All TV values are, hence, normalized to the value obtained with the FIT on the finest, nonadaptive grid, which has been employed. For comparability, the settings were chosen such that all simulations finish within less than one hour.
A comparison of the results shows that: -The relative errors of the FIT and the LT-FIT are similar, however, the TV of the LT-FIT solution is significantly lower, indicating reduced oscillations due to better numerical dispersion properties.
-For identical numbers of DoF the LT-FIT simulations take longer because of the higher computational costs for the time integration.
-For identical mesh resolutions in the bunch area the adaptive simulations yield errors similar to those obtained on nonadaptive meshes. This follows from a comparison of the result of the adaptive simulations with the respective nonadaptive simulation.
Hence, the accuracy of the solutions for this example is mainly determined by the resolution in the bunch region.

B. Self-consistent simulation of the PITZ RF gun
In order to drive a free-electron-laser (FEL) operating by the self-amplified spontaneousemission (SASE) principle, highly charged electron beams of high brightness are required [23].
The aim of the PITZ project (Photo Injector Test Facility at DESY Zeuthen) [24] is the development and testing of an injector capable of delivering such high quality beams. The injectors for the Free Electron Laser in Hamburg (FLASH) [25] and the future European X-Ray Laser Project XFEL [26] are under development at PITZ.
The layout of the radio-frequency (RF) gun of the injector is shown in Fig. 11. The  Table IV.
Many numerical studies of the injector using PIC codes have been performed over the last years. Results have been published, e.g., in [28][29][30][31][32]. However, these simulations either assumed a rotationally symmetric geometry and made use of a so-called 2.5-dimensional approach as implemented in the MAFIA TS2 module [33], or otherwise very short distances in the cm-range of the full three-dimensional model were simulated. Results of parallelized simulations of the full three-dimensional model up to one meter downstream of the cathode have been presented in [4,34]. However, the exact position of the transverse emittance minimum is required for the optimal placement of the accelerating cavity. Also, the value of the emittance at this position is of strong interest in order to check whether the design value is met.
In the context of this work, the gun was simulated up to a distance of two meters downstream of the cathode. First, settings of the simulation parameter such as the grid resolution and the number of macro particles have been determined. Then, a design study was performed, which addresses the effects of individual injector elements on the beam quality.
Previous investigations have shown that the length of the simulated bunches critically depends on the longitudinal grid step size [35]. In [28] it was stated that a grid resolution of 20 µm longitudinally is required in order to obtain accurate results. In the Fig. 12 After emission the electrons have a very low energy of ≈ 5 eV and space charge forces have a strong influence. Since the particles gain energy quickly, it is sufficient to apply the smallest grid step size only in the immediate vicinity of the cathode. Hence, an additional static grid refinement is applied in this area. In the Fig. 12(b) results for this approach are shown.
A grid with a uniform step size of 80 µm is statically refined within the first centimeter from the cathode. Refinement levels from one to four have been applied, resulting in step sizes of 40, 20, 10, and 5 µm within the refined region. The results are, except for minor discrepancies, identical to those obtained with an equidistant grid of the same minimum step size. The results given in the following were obtained with the LT-FIT method. In Fig. 13, the evolution of the horizontal root-mean-square (RMS) bunch width σ x (a), the RMS bunch length σ z (b), the average particle energy (c), and the horizontal emittance ε x (d) for the model shown in Fig. 11 are plotted. The projected emittances ε x and ε y are given by ε u = σ 2 u · σ 2 pu − (σ u,pu ) 2 1 /2 , u ∈ {x, y}, where σ 2 pu is the momentum variance and σ u,pu is the covariance of position and momentum.
We performed a design study of the injector in order to identify the individual effects of the diagnostics section, the laser mirror, and the shutter valve on the beam quality in terms of emittance growth. The respective models are shown in the Figs. 14(a)-(d). For each model, an additional element is added. The colors of the curves in Fig. 14 For the model given in Fig. 14(a), the horizontal and vertical emittances, ε x and ε y , are identical. This is expected since the geometry is rotationally symmetric for all parts close to the beam. The RF input coupler is obviously not rotationally symmetric but this symmetry violation is hidden by the coupling antenna. In Fig. 14(b), the large opening on the bottom side of the doublecross introduces an asymmetry of the geometry in horizontal and vertical direction. While ε x is almost unchanged, the value of ε y increases by approximately 0.04 mm mrad. The laser mirror, added in the horizontal plane, introduces another asymmetry in Fig. 14(c). An influence on the transverse emittances can be identified in the results but the actual emittance growth at the position of the minimum is small. Finally, the shutter valve depicted in Fig. 14(d) results in significant changes of the transverse emittances. While the horizontal emittance actually decreases by approximately 0.1 mm mrad, the vertical emittance increases by about 0.18 mm mrad in comparison to the structure of Fig. 14(a). Hence, the shutter valve causes a distinct asymmetry in the transverse beam dynamics. In the meantime, a shutter valve with RF shieldings has been installed.

VII. CONCLUSIONS
In this paper, we extended the framework of Finite Integration Technique to include dynamic mesh refinement. This offers computability on workstations for a class of multiscale problems, which before was accessible to massively parallelized computations only. In particular, our approach enabled us to perform fully self-consistent simulations of the PITZ RF gun up to two meters downstream from the cathode. We elaborated on the details concerning the mesh refinement and coarsening procedures and presented a novel type of sub-spline interpolation, which is entirely devoid of overshooting. Finally, we performed a design study, which identifies the emittance growth due to individual parts of the gun.

ACKNOWLEDGMENTS
The work of S. Schnepp is supported by the 'Initiative for Excellence' of the German